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

[Satoshi Kimura]

Adding iceshelf parameterization to the trunk.
As of today, this branch has green light from Buildbot (http://buildbot-ocean.ese.ic.ac.uk:8080/builders/iceshelf/builds/9).
The merge will enable us to:
1) prescribe ice as a boundary condition on T and S;
2) allow inhomogeneous Neumann BC on scalar fields discretized in P1DG.

[Alexandros Avdis (alexandros-avdis)]

The in-homogeneous Neumann BC is incomplete, will be making some additions (really made by Stefan) to lp:~fluidity-core/fluidity/DG_Neumann.

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

4021. By Jon Hill on 2012-08-02

New GLS updates which include:
 - Ocean parameterisation which parameterises internal wave breaking, adding some surface TKE into deeper water
 - Updates to assembly routines which are now mostly element-wise
 - Updates to tests
 - Manual updates

4022. By Jon Hill on 2012-08-02

Adding tests for PAR and updating manual

4023. By Rhodri Davies on 2012-08-06

Bug fix for CV, when a diffusion term is present. Again supplied by Cian.

4024. By Christian Jacobs on 2012-08-09

Merging compressible-multiphase branch revisions r3955 to r4067 (inclusive) into trunk.

Summary of changes:
===================

- Added multiphase support for the compressible projection method (currently only for CG Pressure).

- Small changes to the logic in Momentum_Equation.F90. Replaced the use_compressible_projection logical with compressible_eos which is .true. if option_count("/material_phase/equation_of_state/compressible") > 0, and .false. otherwise.

- Added multiphase support for the InternalEnergy equation. Note that the whole state array is now passed to solve_field_equation_cg.

- Implemented the heat transfer term by Gunn (1978).

- Added three P2-P1 multiphase MMS tests for model verification: one without InternalEnergy fields, one with InternalEnergy fields, and one with InternalEnergy fields and heat transfer. The tests will become longtests after merging into the trunk.

- The manual now has documentation on the compressible multiphase model, including the multiphase version of the InternalEnergy equation.

- Added a diagnostic field called CompressibleContinuityResidual which computes the residual of the compressible multiphase continuity equation. This is similar to the SumVelocityDivergence field for incompressible multiphase flow which computes \sum{div(vfrac*u)}.

- Added two gas-solid shock tube tests to help to validate the model: mphase_rogue_shock_tube_dense_bed_glass and mphase_rogue_shock_tube_dense_bed_nylon.

- Added another gas-solid shock tube test (mphase_dusty_gas_shock_tube) which uses the setup by Miura & Glass (1982).

- Implemented two new correlations for fluid-particle drag. The one by Wen & Yu (1966) can support higher Reynolds numbers (compared to the existing Stokes' law correlation). The one by Ergun (1952) can support dense flows where the volume fraction of the particle phase is greater than 0.2, like in the dense bed shock tube tests (and in the pyroclastic flow tests which will hopefully come later).

- Changed the structure of the schema. Users now explicitly select which interaction terms they wish to include, as well as the drag correlation, under "/multiphase_interaction" in Diamond. These options are stored in the new files multiphase_interaction.rnc and multiphase_interaction.rng.

- Added a test for the Wen & Yu (1966) drag correlation. This is similar to mphase_stokes_law.

- Added a test to check if inlet velocity boundary conditions work with the incompressible multiphase flow model.

- Added a test for the flux boundary condition.

4025. By Christian Jacobs on 2012-08-09

Moved the compressible multiphase MMS tests to the longtests directory.

4026. By Christian Jacobs on 2012-08-09

- Updated the tephra_settling example's .flml file to use the new multiphase_interaction schema structure.
- Corrected a small typo in the manual.

4027. By Lawrence Mitchell on 2012-08-13

Fix DG momentum assembly corner case

Since r3984 we would assemble elements in DG momentum if
element_neighbour_owned returned true. However, for some corner cases
element_neighbour_owned need not return true if the element is owned
(consider the case of a process with one owned element). In this
situation we would not correctly assemble all elements. To fix this,
change the test for element assembly to directly check for element
ownership.

At the same time, update the documentation for element_neighbour_owned
to note this issue.

4028. By Rhodri Davies on 2012-08-16

Merge relax_halo_verifies branch into trunk. Thanks to dham for review.

4029. By j-bull08 <email address hidden> on 2012-08-17

Removing duplicate file in 3D BFS example.

4030. By Lawrence Mitchell on 2012-08-17

Fix CG momentum assembly on intel fortran

This only affects the have_les and .not. dynamic_les code path.

4031. By Rhodri Davies on 2012-08-17

Merge zoltan configure check for .mod files into trunk. Thanks to steph for help and jon for review.

4032. By Christian Jacobs on 2012-08-17

Added solvers_converged pass_tests to all multiphase tests.

4033. By Christian Jacobs on 2012-08-19

Added a solvers_converged pass_test to the tephra_settling example.

4034. By Samuel Parkinson on 2012-08-19

Updated Momentum_CG.F90 so that upwind stabilisation can be used with partial or full stress form. The issue is that when using partial or full stress form we tend to set all components of viscosity to the same value. This makes the tensor determinant zero and we cannot invert it. This fix simply sets all non-diagonal terms to zero for the purpose of handling stabilisation when using partial or full stress forms. This isn't ideal as it ignores the fact that we may have non-isotropic viscosities, but these aren't handled properly for these stress forms anyway - so this is a bigger issue that may need fixing later.

The same issue arises when we try to calculate the GridReynoldsNumber diagnostic field and the same fix has been applied here.

4035. By Samuel Parkinson on 2012-08-21

bug fix of SUPG stabilisation in Momentum_CG.F90

[Tim Greaves (tim-greaves)]

Toshi is keen to get reviews of this and make any corrections needed prior to merging; if anyone with knowledge of this code is able to look through this request and give feedback that would be much appreciated!

[Jon Hill (jon-hill)]

Also note that Toshi has made a change to the VTK interface which from the code looks like it allows the user to stop the DG-ifiying of output fields. Are the DG fields projected to continuous somewhere?

Are the fluidity.log and .err files meant to be there?

Looks OK to me, but I would prefer someone else cast an eye over it too.

[Stephan Kramer (s-kramer)]

The DG Neumann bc is being further developed on this branch https://code.launchpad.net/~fluidity-core/fluidity/DG_Neumann. It needed fixing to work with subcycling. It now needs some proper tests.

This merge seems to have all sorts of other unrelated stuff: FLHere(), /io/adhere_to_output_mesh_continuity, /top_position/align_with_geoids, regularise_aspect_ratio, subshelf_hydrostatic_pressure... These all may, or may not be useful additions but should be proposed for a merge seperately.

Please have a look at the diff with the trunk yourself before you put anything up for merge proposal. Reviewing is a lot of work, so we expect things to be cleaned up properly before requesting a merge (no commented out code, etc.)

review: Disapprove

[Adam Candy (asc)]

The latest commit deleting the ice shelf test cases to get the build machine passing was a retrograde step! These need to be fixed and returned to the merge.

review: Needs Fixing

[Adam Candy (asc)]

This branch also contains features currently being worked on in other branches, which themselves are under merge review. These features are not required by the ice shelf thermodynamics, which is what this merge concerns and should be removed from the branch.

These include:
lp:~asc/fluidity/extrudeup
lp:~asc/fluidity/regularisespatialaspect
lp:~asc/fluidity/here
lp:~asc/fluidity/dgifychoice

review: Needs Fixing

[Satoshi Kimura (skimura04)]

Thanks everyone for good thoughts.
I will start from merging in the thermodynamics part first and work the way up.

Toshi

> This branch also contains features currently being worked on in other
> branches, which themselves are under merge review. These features are not
> required by the ice shelf thermodynamics, which is what this merge concerns
> and should be removed from the branch.
>
> These include:
> lp:~asc/fluidity/extrudeup
> lp:~asc/fluidity/regularisespatialaspect
> lp:~asc/fluidity/here
> lp:~asc/fluidity/dgifychoice

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

4036. By Tim Greaves on 2012-08-22

One-line fix, setting the correct flml filename in the test xml.

4037. By Adam Candy on 2012-08-24

Removal of a runtime-generated test file.
jhill - please check!

4038. By Samuel Parkinson on 2012-08-28

Added options for how the bed shear stress is calculated. Yuo can now either:
- use drag formulation as before
- calculate the shear stress based upon the gradient of the velocity field at the bed
- use the shear stress estimate generated by the wall boundary layer boundary condition algorithm

4039. By Samuel Parkinson on 2012-08-31

change to lock_exchange example so that, when run as a test, it runs to 30 seconds so that the tests are checked properly

4040. By Christian Jacobs on 2012-08-31

Merging compressible-multiphase branch revisions r4068 to r4077 (inclusive) into trunk.

Summary of changes:
===================

- Added multiphase support for stress_form, partial_stress_form and isotropic diagonal_viscosity in Momentum_CG.F90. The stress_form will be important in compressible multiphase flow simulations.

- Added an MMS test (mphase_mms_p2p1_stress_form) for the stress_form changes.

4041. By Christian Jacobs on 2012-08-31

Small edit to the manual to generalise the form of the multiphase stress term.

4042. By Xiaohu Guo on 2012-09-01

  1. thread matrix assembly part for Advection_Diffusion_DG.F90
  2. MASSLUMPED_RT0 diffusion scheme is not supported, simply set threads back to 1 if MASSLUMPED_RT0 diffusion scheme is being used, and there will be a warning message.
  3. thread matrix assembly part for Advection_Diffusion_FV.F90

4043. By Christian Jacobs on 2012-09-02

Merging k_epsilon_extend_and_fix branch into trunk.

Changes to k-epsilon model:
===========================
- The code no longer lumps the mass matrix when calculating source terms for the k and epsilon fields, allowing element pairs other than P1-P1 to be used.
- k-epsilon field calculation moved to the beginning of the time-step and calculated with all other diagnostic fields.
- Added source term to momentum equation (-2/3 grad(rho*k)).
- Updated the manual documentation for the k-epsilon model.
- Users must now use partial-stress or full-stress form for the stress term in the momentum equation. This is the correct formulation when using spatially
varying viscosities. The ScalarEddyViscosity from the model is added to all components of the stress tensor to reflect this.
- Fixed calculation of c_eps_3.
- Added lots of debugging options.
- Improved the options_check subroutine.
- Changed how prescribed source terms are implemented for the k and epsilon fields so that multiple non-linear iterations can be carried out
with prescribed source terms.
- Altered calculation of Reynolds Stresses.
- Added a new equation_type (KEpsilon) for the k-epsilon equation. This includes the density in the advection-diffusion equation for the k and epsilon fields. Also made many other changes to the code to take account of density variations. Supports CG and CV discretisations.
- Made a few preliminary steps for adding multiphase support to the k-epsilon model.
- k and epsilon fields are now clipped at Gauss points rather than nodes.
- Moved element assembly to its own subroutine to avoid dynamic memory assignment.
- Changed k and epsilon diffusivity tensor to (nu + nu_t/sigma_*) instead of (nu_t/sigma_*)
- The turbulent buoyancy term is now calculated from the VelocityBuoyancyDensity rather than a summation of scalar fields that have been identified as
being buoyant. This means that any buoyant field will automatically result in a buoyancy term in the k-epsilon model (unless the entire term
is disabled in debugging_options)
- Removed several unused variables.
- Added the buoyancy term to the k-epsilon model.
- Added damping functions for the low_Re k-epsilon model.
- Changed epsilon low_Re boundary condition to d(eps)/dx = 0 as it's simpler than what we had before and still valid (Wilcox 1998).
- Fixes to high-Re boundary conditions.
- Removed length scale limiting.

Changes to tests and examples:
==============================
- The existing k-epsilon MMS test (k-epsilon_mms_cg_p1p1) has been removed and replaced with one that uses a P2-P1 discretisation and runs to steady state.
- Three additional tests have been added for the low_Re boundary condition, and the LinearMomentum Velocity equation_type with CG and CV discretisations.
- Tests include the turbulent buoyancy term and scalar field turbulent diffusion.
- A unit test has been added for the tensor_inner_product function.
- Fixed the backward facing step 2D example based on the changes to the k-epsilon model.

4044. By Christian Jacobs on 2012-09-02

Moved tests/mms_rans_p2p1_keps_linearmomentum and tests/mms_rans_p2p1_keps_linearmomentum_cv to the longtests directory.

4045. By Christian Jacobs on 2012-09-03

Buildbot's makefiles-x86_64 queue has been failing since r4042. These updates from make makefiles will hopefully fix this.

4046. By Adam Candy on 2012-09-03

Wetting and Drying bug fix - the material phase had an assumed name of 'water'.
(It is now just assumed to be the first material phase)

4047. By Satoshi Kimura on 2012-09-04

Added the meltrate code

4048. By Satoshi Kimura on 2012-09-04

Fixed the makefile.dependencies

Unmerged revisions

4058. By Satoshi Kimura on 2013-10-25

Merge the trunk

4057. By Satoshi Kimura on 2013-06-11

Merge the trunk

4056. By Satoshi Kimura on 2013-02-20

Merge the recent trunk

4055. By Satoshi Kimura on 2013-02-19

Changed the min velocity for melting

4054. By Satoshi Kimura on 2013-02-19

Changed the min velocity for melting

4053. By Satoshi Kimura on 2012-10-30

Merged the trunk to be updated

4052. By Satoshi Kimura on 2012-09-06

Fixing the test case, iceshelf

4051. By Satoshi Kimura on 2012-09-06

Fixed my test case to pass the medium test

4050. By Satoshi Kimura on 2012-09-05

Deleted some comments in Ice_Melt_Interface.F90

4049. By Satoshi Kimura on 2012-09-05

Added the test case

=== modified file 'aclocal.m4'
--- aclocal.m4	2012-08-16 17:13:12 +0000
+++ aclocal.m4	2012-09-04 15:00:37 +0000
@@ -412,6 +412,8 @@
   AC_MSG_ERROR( [You need to set PETSC_DIR to point at your PETSc installation... exiting] )
 fi
 
+
+
 PETSC_LINK_LIBS=`make -s -f petsc_makefile getlinklibs`
 LIBS="$PETSC_LINK_LIBS $LIBS"
 
@@ -465,7 +467,7 @@
 # petsc tutorial - using the headers in the same way as we do in the code
 AC_LINK_IFELSE(
 [AC_LANG_SOURCE([
-program test_petsc
+program test
 #include "petscversion.h"
 #ifdef HAVE_PETSC_MODULES
   use petsc
@@ -572,13 +574,13 @@
       call MatDestroy(A,ierr)
 
       call PetscFinalize(ierr)
-end program test_petsc
+end program test
 ])],
 [
 AC_MSG_NOTICE([PETSc program succesfully compiled and linked.])
 ],
 [
-cp conftest.F90 test_petsc.F90
+cp conftest.F90 test.F90
 AC_MSG_FAILURE([Failed to compile and link PETSc program.])])
 
 AC_LANG_RESTORE
@@ -1004,6 +1006,7 @@
 	echo "note: using $zoltan_INCLUDES_PATH/zoltan.mod"
 fi 
 
+
 AC_LANG_SAVE
 AC_LANG_C
 AC_CHECK_LIB(
@@ -1011,28 +1014,7 @@
 	[Zoltan_Initialize],
 	[AC_DEFINE(HAVE_ZOLTAN,1,[Define if you have zoltan library.])],
 	[AC_MSG_ERROR( [Could not link in the zoltan library... exiting] )] )
-
-# Small test for zoltan .mod files:
-AC_LANG(Fortran)
-ac_ext=F90
-# In fluidity's makefile we explicitly add CPPFLAGS, temporarily add it to
-# FCFLAGS here for this zoltan test:
-tmpFCFLAGS="$FCFLAGS"
-FCFLAGS="$FCFLAGS $CPPFLAGS"
-AC_LINK_IFELSE(
-[AC_LANG_SOURCE([
-program test_zoltan
- use zoltan
-end program test_zoltan
-])],
-[
-AC_MSG_NOTICE([Great success! Zoltan .mod files exist and are usable])
-],
-[
-cp conftest.F90 test_zoltan.F90
-AC_MSG_FAILURE([Failed to find zoltan.mod files])])
-# And now revert FCFLAGS
-FCFLAGS="$tmpFCFLAGS"
+# Save variables...
 AC_LANG_RESTORE
 
 ZOLTAN="yes"

=== modified file 'configure'
--- configure	2012-08-16 17:13:12 +0000
+++ configure	2012-09-04 15:00:37 +0000
@@ -2137,52 +2137,6 @@
 
 } # ac_fn_f77_try_link
 
-# ac_fn_fc_try_link LINENO
-# ------------------------
-# Try to link conftest.$ac_ext, and return whether this succeeded.
-ac_fn_fc_try_link ()
-{
-  as_lineno=${as_lineno-"$1"} as_lineno_stack=as_lineno_stack=$as_lineno_stack
-  rm -f conftest.$ac_objext conftest$ac_exeext
-  if { { ac_try="$ac_link"
-case "(($ac_try" in
-  *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;;
-  *) ac_try_echo=$ac_try;;
-esac
-eval ac_try_echo="\"\$as_me:${as_lineno-$LINENO}: $ac_try_echo\""
-$as_echo "$ac_try_echo"; } >&5
-  (eval "$ac_link") 2>conftest.err
-  ac_status=$?
-  if test -s conftest.err; then
-    grep -v '^ *+' conftest.err >conftest.er1
-    cat conftest.er1 >&5
-    mv -f conftest.er1 conftest.err
-  fi
-  $as_echo "$as_me:${as_lineno-$LINENO}: \$? = $ac_status" >&5
-  test $ac_status = 0; } && {
-	 test -z "$ac_fc_werror_flag" ||
-	 test ! -s conftest.err
-       } && test -s conftest$ac_exeext && {
-	 test "$cross_compiling" = yes ||
-	 $as_test_x conftest$ac_exeext
-       }; then :
-  ac_retval=0
-else
-  $as_echo "$as_me: failed program was:" >&5
-sed 's/^/| /' conftest.$ac_ext >&5
-
-	ac_retval=1
-fi
-  # Delete the IPA/IPO (Inter Procedural Analysis/Optimization) information
-  # created by the PGI compiler (conftest_ipa8_conftest.oo), as it would
-  # interfere with the next link command; also delete a directory that is
-  # left behind by Apple's compiler.  We do this before executing the actions.
-  rm -rf conftest.dSYM conftest_ipa8_conftest.oo
-  eval $as_lineno_stack; ${as_lineno_stack:+:} unset as_lineno
-  as_fn_set_status $ac_retval
-
-} # ac_fn_fc_try_link
-
 # ac_fn_c_check_func LINENO FUNC VAR
 # ----------------------------------
 # Tests whether FUNC exists, setting the cache variable VAR accordingly
@@ -2266,6 +2220,52 @@
 
 } # ac_fn_c_check_func
 
+# ac_fn_fc_try_link LINENO
+# ------------------------
+# Try to link conftest.$ac_ext, and return whether this succeeded.
+ac_fn_fc_try_link ()
+{
+  as_lineno=${as_lineno-"$1"} as_lineno_stack=as_lineno_stack=$as_lineno_stack
+  rm -f conftest.$ac_objext conftest$ac_exeext
+  if { { ac_try="$ac_link"
+case "(($ac_try" in
+  *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;;
+  *) ac_try_echo=$ac_try;;
+esac
+eval ac_try_echo="\"\$as_me:${as_lineno-$LINENO}: $ac_try_echo\""
+$as_echo "$ac_try_echo"; } >&5
+  (eval "$ac_link") 2>conftest.err
+  ac_status=$?
+  if test -s conftest.err; then
+    grep -v '^ *+' conftest.err >conftest.er1
+    cat conftest.er1 >&5
+    mv -f conftest.er1 conftest.err
+  fi
+  $as_echo "$as_me:${as_lineno-$LINENO}: \$? = $ac_status" >&5
+  test $ac_status = 0; } && {
+	 test -z "$ac_fc_werror_flag" ||
+	 test ! -s conftest.err
+       } && test -s conftest$ac_exeext && {
+	 test "$cross_compiling" = yes ||
+	 $as_test_x conftest$ac_exeext
+       }; then :
+  ac_retval=0
+else
+  $as_echo "$as_me: failed program was:" >&5
+sed 's/^/| /' conftest.$ac_ext >&5
+
+	ac_retval=1
+fi
+  # Delete the IPA/IPO (Inter Procedural Analysis/Optimization) information
+  # created by the PGI compiler (conftest_ipa8_conftest.oo), as it would
+  # interfere with the next link command; also delete a directory that is
+  # left behind by Apple's compiler.  We do this before executing the actions.
+  rm -rf conftest.dSYM conftest_ipa8_conftest.oo
+  eval $as_lineno_stack; ${as_lineno_stack:+:} unset as_lineno
+  as_fn_set_status $ac_retval
+
+} # ac_fn_fc_try_link
+
 # ac_fn_fc_try_run LINENO
 # -----------------------
 # Try to link conftest.$ac_ext, and return whether this succeeded. Assumes
@@ -11791,6 +11791,7 @@
 fi
 
 
+
 ac_ext=c
 ac_cpp='$CPP $CPPFLAGS'
 ac_compile='$CC -c $CFLAGS $CPPFLAGS conftest.$ac_ext >&5'
@@ -11857,42 +11858,7 @@
   as_fn_error $? "Could not link in the zoltan library... exiting " "$LINENO" 5
 fi
 
-
-# Small test for zoltan .mod files:
-ac_ext=${ac_fc_srcext-f}
-ac_compile='$FC -c $FCFLAGS $ac_fcflags_srcext conftest.$ac_ext >&5'
-ac_link='$FC -o conftest$ac_exeext $FCFLAGS $LDFLAGS $ac_fcflags_srcext conftest.$ac_ext $LIBS >&5'
-ac_compiler_gnu=$ac_cv_fc_compiler_gnu
-
-ac_ext=F90
-# In fluidity's makefile we explicitly add CPPFLAGS, temporarily add it to
-# FCFLAGS here for this zoltan test:
-tmpFCFLAGS="$FCFLAGS"
-FCFLAGS="$FCFLAGS $CPPFLAGS"
-cat > conftest.$ac_ext <<_ACEOF
-
-program test_zoltan
- use zoltan
-end program test_zoltan
-
-_ACEOF
-if ac_fn_fc_try_link "$LINENO"; then :
-
-{ $as_echo "$as_me:${as_lineno-$LINENO}: Great success! Zoltan .mod files exist and are usable" >&5
-$as_echo "$as_me: Great success! Zoltan .mod files exist and are usable" >&6;}
-
-else
-
-cp conftest.F90 test_zoltan.F90
-{ { $as_echo "$as_me:${as_lineno-$LINENO}: error: in \`$ac_pwd':" >&5
-$as_echo "$as_me: error: in \`$ac_pwd':" >&2;}
-as_fn_error $? "Failed to find zoltan.mod files
-See \`config.log' for more details" "$LINENO" 5; }
-fi
-rm -f core conftest.err conftest.$ac_objext \
-    conftest$ac_exeext conftest.$ac_ext
-# And now revert FCFLAGS
-FCFLAGS="$tmpFCFLAGS"
+# Save variables...
 ac_ext=c
 ac_cpp='$CPP $CPPFLAGS'
 ac_compile='$CC -c $CFLAGS $CPPFLAGS conftest.$ac_ext >&5'
@@ -12935,6 +12901,8 @@
   as_fn_error $? "You need to set PETSC_DIR to point at your PETSc installation... exiting " "$LINENO" 5
 fi
 
+
+
 PETSC_LINK_LIBS=`make -s -f petsc_makefile getlinklibs`
 LIBS="$PETSC_LINK_LIBS $LIBS"
 
@@ -13004,7 +12972,7 @@
 # petsc tutorial - using the headers in the same way as we do in the code
 cat > conftest.$ac_ext <<_ACEOF
 
-program test_petsc
+program test
 #include "petscversion.h"
 #ifdef HAVE_PETSC_MODULES
   use petsc
@@ -13111,7 +13079,7 @@
       call MatDestroy(A,ierr)
 
       call PetscFinalize(ierr)
-end program test_petsc
+end program test
 
 _ACEOF
 if ac_fn_fc_try_link "$LINENO"; then :
@@ -13121,7 +13089,7 @@
 
 else
 
-cp conftest.F90 test_petsc.F90
+cp conftest.F90 test.F90
 { { $as_echo "$as_me:${as_lineno-$LINENO}: error: in \`$ac_pwd':" >&5
 $as_echo "$as_me: error: in \`$ac_pwd':" >&2;}
 as_fn_error $? "Failed to compile and link PETSc program.

=== added file 'examples/backward_facing_step_3d/Ercoftac-test31-BFS/readme'
--- examples/backward_facing_step_3d/Ercoftac-test31-BFS/readme	1970-01-01 00:00:00 +0000
+++ examples/backward_facing_step_3d/Ercoftac-test31-BFS/readme	2012-09-04 15:00:37 +0000
@@ -0,0 +1,28 @@
+1) File name convention:  filename.nnn
+   nnn = streamwise nodal (cell center) location
+  
+   nnn = 181  --->  x/h = -3.0
+   nnn = 360  --->  x/h = 4.0
+   nnn = 411  --->  x/h = 6.0
+   nnn = 513  --->  x/h = 10.0
+   nnn = 641  --->  x/h = 15.0
+   nnn = 744  --->  x/h = 19.0
+
+2) Files "x.nnn" contain U, V, u'(rms), v'(rms), w'(rms), u'v'.  
+   File header explains the data columns.
+
+   All quantities are nomalized to U0, where U0 is the mean inlet 
+   free stream velocity.
+
+3) Files "r**.nnn" contain Reynolds stress budgets.
+
+   r11.nnn  --->  Streamwise component (u'u')
+   r22.nnn  --->  Vertical component (v'v')
+   r33.nnn  --->  Spanwise component (w'w')
+   r12.nnn  --->  Reynolds shear stress component (u'v')
+   rkk.nnn  --->  Turbulent kinetic energy
+
+   All quantities are nomalized to U0^3/h.
+
+4) File "stat.info" contains additional information at each streamwise 
+   location (u_tau, Cf, etc...).

=== modified file 'examples/backward_facing_step_3d/Makefile'
--- examples/backward_facing_step_3d/Makefile	2012-07-27 08:10:27 +0000
+++ examples/backward_facing_step_3d/Makefile	2012-09-04 15:00:37 +0000
@@ -39,7 +39,7 @@
 	./postprocessor_3d.py
 
 input: clean
-	$(MAKE) preprocess NPROCS=8
+	$(MAKE) preprocess
 
 clean:
 	@echo **********Cleaning the output from previous fluidity runs:

=== modified file 'examples/hokkaido-nansei-oki_tsunami/hokkaido-nansei-oki_tsunami.xml'
--- examples/hokkaido-nansei-oki_tsunami/hokkaido-nansei-oki_tsunami.xml	2012-08-22 21:21:08 +0000
+++ examples/hokkaido-nansei-oki_tsunami/hokkaido-nansei-oki_tsunami.xml	2012-09-04 15:00:37 +0000
@@ -2,8 +2,8 @@
 <testproblem>
   <name>Monai Valley</name>
   <owner userid="sf1409"/>
-  <problem_definition length="long" nprocs="4">
-    <command_line>fluidity MonaiValley_C_p1p1_nu0.01_kmkstab_drag0.002_butcircularoundisland0.2.flml</command_line>
+  <problem_definition length="special" nprocs="4">
+    <command_line>fluidity MonaiValley_C.flml</command_line>
   </problem_definition>
   <variables>
     <variable name="gage_integral_error" language="python">import tools

=== modified file 'examples/restratification_after_oodc/Makefile'
--- examples/restratification_after_oodc/Makefile	2012-07-25 15:27:53 +0000
+++ examples/restratification_after_oodc/Makefile	2012-09-04 15:00:37 +0000
@@ -40,7 +40,7 @@
 	@echo **********No postprocessing needed
 
 input: clean
-	$(MAKE) preprocess NPROCS=32
+	$(MAKE) preprocess NPROCS=8
 
 clean:
 	@echo **********Cleaning the output from previous fluidity runs:

=== modified file 'examples/tephra_settling/tephra_settling.flml'
--- examples/tephra_settling/tephra_settling.flml	2012-08-09 22:08:12 +0000
+++ examples/tephra_settling/tephra_settling.flml	2012-09-04 15:00:37 +0000
@@ -586,9 +586,4 @@
       </particle_diameter>
     </multiphase_properties>
   </material_phase>
-  <multiphase_interaction>
-    <fluid_particle_drag>
-      <drag_correlation name="stokes"/>
-    </fluid_particle_drag>
-  </multiphase_interaction>
 </fluidity_options>

=== modified file 'examples/tephra_settling/tephra_settling.xml'
--- examples/tephra_settling/tephra_settling.xml	2012-08-19 00:21:26 +0000
+++ examples/tephra_settling/tephra_settling.xml	2012-09-04 15:00:37 +0000
@@ -20,18 +20,11 @@
 s = stat_parser("tephra_settling.stat")
 tephra_u_max=s["Tephra"]["Velocity%magnitude"]["max"]
     </variable>
-    
     <variable name="time" language="python">
 from fluidity_tools import stat_parser
 s = stat_parser("tephra_settling.stat")
 time=s["ElapsedTime"]["value"]
     </variable>
-    
-    <variable name="solvers_converged" language="python">
-import os
-files = os.listdir("./")
-solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files
-    </variable>
   </variables>
     
   <pass_tests>
@@ -49,10 +42,6 @@
       assert max(tephra_u_max[t:]) &lt; 0.05
       break
     </test>
-    
-    <test name="Solvers converged" language="python">
-assert(solvers_converged)
-    </test>
   </pass_tests>
   
   <warn_tests>

=== modified file 'examples/top_hat/top_hat.xml'
--- examples/top_hat/top_hat.xml	2012-07-25 12:38:28 +0000
+++ examples/top_hat/top_hat.xml	2012-09-04 15:00:37 +0000
@@ -2,7 +2,7 @@
 <testproblem>
   <name>tophat</name>
   <owner userid="skramer"/>
-  <problem_definition length="long" nprocs="1">
+  <problem_definition length="medium" nprocs="1">
     <command_line>fluidity -v2 -l top_hat_cg.flml &amp;&amp; 
 fluidity -v2 -l top_hat_dg.flml &amp;&amp;
 fluidity -v2 -l top_hat_cv.flml</command_line>

=== modified file 'femtools/Diagnostic_Fields.F90'
--- femtools/Diagnostic_Fields.F90	2012-08-29 18:09:06 +0000
+++ femtools/Diagnostic_Fields.F90	2012-09-04 15:00:37 +0000
@@ -477,7 +477,7 @@
     type(scalar_field), intent(inout) :: grn
 
     type(vector_field), pointer :: U, X
-    integer :: ele, gi, stat, a, b
+    integer :: ele, gi, stat
     ! Transformed quadrature weights.
     real, dimension(ele_ngi(GRN, 1)) :: detwei
     ! Inverse of the local coordinate change matrix.
@@ -525,24 +525,6 @@
        ! The matmul is as given by dham
        GRN_q=ele_val_at_quad(U, ele)
        vis_q=ele_val_at_quad(viscosity, ele)
-
-       ! for full and partial stress form we need to set the off diagonal terms of the viscosity tensor to zero
-       ! to be able to invert it 
-       if (have_option(trim(U%option_path)//&
-            &"/prognostic/spatial_discretisation/continuous_galerkin"//&
-            &"/stress_terms/stress_form") .or. &
-            have_option(trim(U%option_path)//&
-            &"/prognostic/spatial_discretisation/continuous_galerkin"//&
-            &"/stress_terms/partial_stress_form")) then
-
-          do a=1,size(vis_q,1)
-             do b=1,size(vis_q,2)
-                if(a.eq.b) cycle
-                vis_q(a,b,:) = 0.0
-             end do
-          end do
-       end if
-
        do gi=1, size(detwei)
           GRN_q(:,gi)=matmul(GRN_q(:,gi), J(:,:,gi))
           GRN_q(:,gi)=matmul(inverse(vis_q(:,:,gi)), GRN_q(:,gi))

=== modified file 'femtools/Halos_Communications.F90'
--- femtools/Halos_Communications.F90	2012-08-10 16:03:41 +0000
+++ femtools/Halos_Communications.F90	2012-09-04 15:00:37 +0000
@@ -765,8 +765,6 @@
     
     type(halo_type), intent(in) :: halo
     real, dimension(:), intent(in) :: real_array
-
-    real :: epsl
     
     logical :: verifies
     
@@ -780,16 +778,14 @@
     
     call halo_update(halo, lreal_array)
     
-    epsl = spacing( maxval( abs( lreal_array ))) * 10000.
-    call allmax(epsl)
-
-    verifies = all(abs(real_array - lreal_array) < epsl )
+    verifies = all(abs(real_array - lreal_array) < 10000.0 * max(spacing(lreal_array), epsilon(lreal_array)))
 #ifdef DDEBUG
     if(.not. verifies) then
       do i = 1, halo_proc_count(halo)
         do j = 1, halo_receive_count(halo, i)
           receive = halo_receive(halo, i, j)
-          if(abs(real_array(receive) - lreal_array(receive)) >= epsl) then
+          if(abs(real_array(receive) - lreal_array(receive)) >= 10000.0&
+               &*max(spacing(real_array(receive)), epsilon(real_array(receive)))) then 
             ewrite(0, *) "Warning: Halo receive ", receive, " for halo " // halo_name(halo) // " failed verification"
             ewrite(0, *) "Reference = ", real_array(receive)
             ewrite(0, *) "Value in verification array = ", lreal_array(receive)
@@ -798,12 +794,13 @@
       end do
       
       do i = 1, size(real_array)
-         if(abs(real_array(i) - lreal_array(i)) >= epsl ) then
-            ewrite(0, *) "Warning: Reference index ", i, " for halo " // halo_name(halo) // " failed verification"
-            ewrite(0, *) "Reference = ", real_array(i)
-            ewrite(0, *) "Value in verification array = ", lreal_array(i)
-         end if
-     end do
+        if(abs(real_array(i) - lreal_array(i)) >= 10000.0&
+               &*max(spacing(real_array(i)), epsilon(real_array(i)))) then
+          ewrite(0, *) "Warning: Reference index ", i, " for halo " // halo_name(halo) // " failed verification"
+          ewrite(0, *) "Reference = ", real_array(i)
+          ewrite(0, *) "Value in verification array = ", lreal_array(i)
+        end if
+      end do
     end if
 #endif
     

=== modified file 'femtools/Makefile.dependencies'
--- femtools/Makefile.dependencies	2012-05-29 13:43:29 +0000
+++ femtools/Makefile.dependencies	2012-09-04 15:00:37 +0000
@@ -1192,18 +1192,6 @@
 
 Vector_set.o ../include/vector_set.mod: Vector_set.F90
 
-../include/vertical_extrapolation_module.mod: Vertical_Extrapolation.o
-	@true
-
-Vertical_Extrapolation.o ../include/vertical_extrapolation_module.mod: \
-   Vertical_Extrapolation.F90 ../include/boundary_conditions.mod \
-   ../include/dynamic_bin_sort_module.mod ../include/elements.mod \
-   ../include/fdebug.h ../include/fields.mod ../include/fldebug.mod \
-   ../include/global_parameters.mod ../include/integer_set_module.mod \
-   ../include/parallel_tools.mod ../include/pickers.mod \
-   ../include/sparse_tools.mod ../include/spud.mod ../include/state_module.mod \
-   ../include/transform_elements.mod ../include/vtk_interfaces.mod
-
 ../include/wandzura_quadrature.mod: Wandzura_Quadrature.o
 	@true
 

=== modified file 'femtools/Makefile.in'
--- femtools/Makefile.in	2012-06-08 20:28:35 +0000
+++ femtools/Makefile.in	2012-09-04 15:00:37 +0000
@@ -104,7 +104,7 @@
   CGAL_Tools_C.o CGAL_Tools.o \
   Rotated_Boundary_Conditions.o MPI_Interfaces.o Parallel_Tools.o \
   Fields_Halos.o Profiler.o Profiler_Fortran.o Streamfunction.o \
-  GMSH_Common.o Read_GMSH.o Write_GMSH.o Mesh_Files.o Vertical_Extrapolation.o
+  GMSH_Common.o Read_GMSH.o Write_GMSH.o Mesh_Files.o
 
 
 # objects to be included in libfemtools:

=== modified file 'femtools/Parallel_fields.F90'
--- femtools/Parallel_fields.F90	2012-08-13 12:39:57 +0000
+++ femtools/Parallel_fields.F90	2012-09-04 15:00:37 +0000
@@ -284,11 +284,8 @@
     !!< Return .true. if ELEMENT_NUMBER has a neighbour in MESH that
     !!< is owned by this process otherwise .false.
 
-    !! Note, you cannot use this function to compute if
-    !! ELEMENT_NUMBER is owned.  Imagine if this is the only owned
-    !! element on a process, then none of the neighbours will be owned.
-    !! You can use this function to compute whether ELEMENT_NUMBER is
-    !! in the L1 element halo.
+    !! Effectively, this computes whether ELEMENT_NUMBER is in the L1
+    !! halo.
     type(mesh_type), intent(in) :: mesh
     integer, intent(in) :: element_number
     logical :: owned
@@ -314,11 +311,8 @@
     !!< Return .true. if ELEMENT_NUMBER has a neighbour in FIELD that
     !!< is owned by this process otherwise .false.
 
-    !! Note, you cannot use this function to compute if
-    !! ELEMENT_NUMBER is owned.  Imagine if this is the only owned
-    !! element on a process, then none of the neighbours will be owned.
-    !! You can use this function to compute whether ELEMENT_NUMBER is
-    !! in the L1 element halo.
+    !! Effectively, this computes whether ELEMENT_NUMBER is in the L1
+    !! halo.
     type(scalar_field), intent(in) :: field
     integer, intent(in) :: element_number
     logical :: owned
@@ -330,11 +324,8 @@
     !!< Return .true. if ELEMENT_NUMBER has a neighbour in FIELD that
     !!< is owned by this process otherwise .false.
 
-    !! Note, you cannot use this function to compute if
-    !! ELEMENT_NUMBER is owned.  Imagine if this is the only owned
-    !! element on a process, then none of the neighbours will be owned.
-    !! You can use this function to compute whether ELEMENT_NUMBER is
-    !! in the L1 element halo.
+    !! Effectively, this computes whether ELEMENT_NUMBER is in the L1
+    !! halo.
     type(vector_field), intent(in) :: field
     integer, intent(in) :: element_number
     logical :: owned
@@ -346,11 +337,8 @@
     !!< Return .true. if ELEMENT_NUMBER has a neighbour in FIELD that
     !!< is owned by this process otherwise .false.
 
-    !! Note, you cannot use this function to compute if
-    !! ELEMENT_NUMBER is owned.  Imagine if this is the only owned
-    !! element on a process, then none of the neighbours will be owned.
-    !! You can use this function to compute whether ELEMENT_NUMBER is
-    !! in the L1 element halo.
+    !! Effectively, this computes whether ELEMENT_NUMBER is in the L1
+    !! halo.
     type(tensor_field), intent(in) :: field
     integer, intent(in) :: element_number
     logical :: owned

=== removed file 'femtools/Vertical_Extrapolation.F90'
--- femtools/Vertical_Extrapolation.F90	2012-05-29 13:43:29 +0000
+++ femtools/Vertical_Extrapolation.F90	1970-01-01 00:00:00 +0000
@@ -1,1290 +0,0 @@
-!    Copyright (C) 2007 Imperial College London and others.
-!
-!    Please see the AUTHORS file in the main source directory for a full list
-!    of copyright holders.
-!
-!    Prof. C Pain
-!    Applied Modelling and Computation Group
-!    Department of Earth Science and Engineering
-!    Imperial College London
-!
-!    [e-mail]
-!
-!    This library is free software; you can redistribute it and/or
-!    modify it under the terms of the GNU Lesser General Public
-!    License as published by the Free Software Foundation,
-!    version 2.1 of the License.
-!
-!    This library is distributed in the hope that it will be useful,
-!    but WITHOUT ANY WARRANTY; without even the implied warranty of
-!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-!    Lesser General Public License for more details.
-!
-!    You should have received a copy of the GNU Lesser General Public
-!    License along with this library; if not, write to the Free Software
-!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
-!    USA
-
-
-#include "confdefs.h"
-#include "fdebug.h"
-
-module vertical_extrapolation_module
-!!< Module containing routines for vertical extrapolation on
-!!< fully unstructured 3D meshes. Also contains routines for updating
-!!< distance to top and bottom fields.
-use fldebug
-use elements
-use sparse_tools
-use fields
-use state_module
-use transform_elements
-use boundary_conditions
-use parallel_tools
-use global_parameters, only: real_4, real_8
-use spud
-use dynamic_bin_sort_module
-use pickers
-use integer_set_module
-use vtk_interfaces
-implicit none
-
-interface VerticalExtrapolation
-  module procedure VerticalExtrapolationScalar, &
-    VerticalExtrapolationMultiple, VerticalExtrapolationVector
-end interface VerticalExtrapolation
-
-interface vertical_integration
-  module procedure vertical_integration_scalar, &
-    vertical_integration_multiple, vertical_integration_vector
-end interface vertical_integration
-
-!! element i is considered to be above j iff the inner product of the face
-!! normal point from i to j and the gravity normal is bigger than this
-real, parameter:: VERTICAL_INTEGRATION_EPS=1.0e-8
-
-! the two hemispheres are both projected to the unit circle, but the projected southern
-! hemisphere will have its origin shifted. The projected hemispheres are slightly bigger
-! than the unit circle, as a bit of overlap with the opposite hemisphere is included around
-! the equator. So this shift needs to be big enough to not make both projected hemispheres
-! overlap
-real, parameter:: HEMISPHERE_SHIFT=10.0
-  
-private
-
-public CalculateTopBottomDistance
-public VerticalExtrapolation, vertical_integration
-public vertical_element_ordering
-public VerticalProlongationOperator
-public vertical_extrapolation_module_check_options
-public compute_face_normal_gravity
-
-contains
-  
-subroutine UpdateDistanceField(state, name, vertical_coordinate)
-  ! This sub calculates the vertical distance to the free surface
-  ! and bottom of the ocean to all nodes. The results are stored
-  ! in the 'DistanceToBottom/FreeSurface' fields from state.
-  type(state_type), intent(inout):: state
-  character(len=*), intent(in):: name
-  type(scalar_field), intent(in):: vertical_coordinate
-
-  ! Local variables
-  type(vector_field), pointer:: positions, vertical_normal
-  type(scalar_field), pointer:: distance
-  
-  integer, pointer, dimension(:):: surface_element_list
-
-  ! the distance field to compute:
-  distance => extract_scalar_field(state, name)
-  ! its first boundary condition is on the related top or bottom mesh
-  call get_boundary_condition(distance, 1, &
-    surface_element_list=surface_element_list)
-  positions => extract_vector_field(state, "Coordinate")
-  vertical_normal => extract_vector_field(state, "GravityDirection")
-  
-  ! in each node of the mesh, set "distance" to the vertical coordinate 
-  ! of this node projected to the above/below surface mesh
-  call VerticalExtrapolation(vertical_coordinate, distance, positions, &
-    vertical_normal, surface_element_list, surface_name=name)
-    
-  ! the distance is then calculated by subtracting its own vertical coordinate
-  call addto(distance, vertical_coordinate, scale=-1.0)
-  
-  if (name=="DistanceToBottom") then
-    ! make distance to bottom positive
-    call scale(distance, -1.0)
-  end if
-  
-end subroutine UpdateDistanceField
-  
-subroutine CalculateTopBottomDistance(state)
-  !! This sub calculates the vertical distance to the free surface
-  !! and bottom of the ocean to all nodes. The results are stored
-  !! in the 'DistanceToBottom/Top' fields from state.
-  type(state_type), intent(inout):: state
-    
-  type(mesh_type), pointer:: xmesh
-  type(scalar_field):: vertical_coordinate
-  
-  xmesh => extract_mesh(state, "CoordinateMesh")
-  call allocate(vertical_coordinate, xmesh, "VerticalCoordinate")
-  call calculate_vertical_coordinate(state, vertical_coordinate)
-  call UpdateDistanceField(state, "DistanceToTop", vertical_coordinate)
-  call UpdateDistanceField(state, "DistanceToBottom", vertical_coordinate)
-  call deallocate(vertical_coordinate)
-
-end subroutine CalculateTopBottomDistance
-  
-subroutine calculate_vertical_coordinate(state, vertical_coordinate)
-  !! Computes a vertical coordinate, i.e. a scalar field such that
-  !! for each 2 nodes above each other, the difference of the field
-  !! in these nodes gives the distance between them.
-  type(state_type), intent(inout):: state
-  type(scalar_field), intent(inout):: vertical_coordinate
-    
-  type(vector_field), pointer:: positions, gravity_normal
-  type(scalar_field):: positions_magnitude
-  
-  positions => extract_vector_field(state, "Coordinate")
-  if(have_option('/geometry/spherical_earth')) then
-    ! use the radius as vertical coordinate
-    ! that is, the l2-norm of the coordinate field
-    positions_magnitude=magnitude(positions)
-    call set(vertical_coordinate, positions_magnitude)
-    call deallocate(positions_magnitude)
-  else
-    gravity_normal => extract_vector_field(state, "GravityDirection")
-    assert(gravity_normal%field_type==FIELD_TYPE_CONSTANT)
-    call inner_product(vertical_coordinate, gravity_normal, positions)
-    ! gravity points down, we want a vertical coordinate that increases upward
-    call scale(vertical_coordinate, -1.0)
-  end if
-  
-end subroutine calculate_vertical_coordinate
-
-subroutine VerticalExtrapolationScalar(from_field, to_field, &
-  positions, vertical_normal, surface_element_list, surface_name)
-  !!< This sub extrapolates the values on a horizontal 2D surface
-  !!< in the vertical direction to 3D fields
-  !! The from_fields should be 3D fields of which only the values on the
-  !! 2D horizontal surface are used.
-  type(scalar_field), intent(in):: from_field
-  !! Resulting extrapolated field. May be the same field or a field on
-  !! a different mesh (different degree).
-  type(scalar_field), intent(inout):: to_field
-  !! positions, and upward normal vector on the whole domain
-  type(vector_field), target, intent(inout):: positions
-  type(vector_field), target, intent(in):: vertical_normal
-  !! the surface elements (faces numbers) that make up the surface
-  integer, dimension(:), intent(in):: surface_element_list
-  !! If provided the projected surface mesh onto horizontal coordinates 
-  !! and its associated rtree/pickers are cached under this name and 
-  !! attached to the 'positions'. In this case when called again with
-  !! the same 'positions' and the same surface_name,
-  !! the same surface_element_list should again be provided.
-  character(len=*), optional, intent(in):: surface_name
-
-  type(scalar_field), dimension(1):: to_fields
-    
-  to_fields=(/ to_field /)
-  call VerticalExtrapolationMultiple( (/ from_field /) , to_fields, &
-      positions, vertical_normal, surface_element_list, &
-      surface_name=surface_name)
-  
-end subroutine VerticalExtrapolationScalar
-
-subroutine VerticalExtrapolationVector(from_field, to_field, &
-  positions, vertical_normal, surface_element_list, surface_name)
-  !!< This sub extrapolates the values on a horizontal 2D surface
-  !!< in the vertical direction to 3D fields
-  !! The from_fields should be 3D fields of which only the values on the
-  !! 2D horizontal surface are used.
-  type(vector_field), intent(in):: from_field
-  !! Resulting extrapolated field. May be the same field or a field on
-  !! a different mesh (different degree).
-  type(vector_field), intent(inout):: to_field
-  !! positions, and upward normal vector on the whole domain
-  type(vector_field), target, intent(inout):: positions
-  type(vector_field), target, intent(in):: vertical_normal
-  !! the surface elements (faces numbers) that make up the surface
-  integer, dimension(:), intent(in):: surface_element_list
-  !! If provided the projected surface mesh onto horizontal coordinates 
-  !! and its associated rtree/pickers are cached under this name and 
-  !! attached to the 'positions'. In this case when called again with
-  !! the same 'positions' and the same surface_name,
-  !! the same surface_element_list should again be provided.
-  character(len=*), optional, intent(in):: surface_name
-  
-  type(scalar_field), dimension(from_field%dim):: from_field_components, to_field_components
-  integer i
-  
-  assert(from_field%dim==to_field%dim)
-  
-  do i=1, from_field%dim
-     from_field_components(i)=extract_scalar_field(from_field, i)
-     to_field_components(i)=extract_scalar_field(to_field, i)
-  end do
-    
-  call VerticalExtrapolationMultiple( from_field_components, to_field_components, &
-        positions, vertical_normal, surface_element_list, &
-        surface_name=surface_name)
-  
-end subroutine VerticalExtrapolationVector
-
-subroutine VerticalExtrapolationMultiple(from_fields, to_fields, &
-  positions, vertical_normal, surface_element_list, surface_name)
-  !!< This sub extrapolates the values on a horizontal 2D surface
-  !!< in the vertical direction to 3D fields
-  !! The from_fields should be 3D fields of which only the values on the
-  !! 2D horizontal surface are used.
-  !! This version takes multiple from_fields at the same time and extrapolates
-  !! to to_fields, such that the surface search only has to be done once. This 
-  !! will only work if all the from_fields are on the same mesh, and on all the 
-  !! to_field are on the same (possibly a different) mesh.
-  type(scalar_field), dimension(:), intent(in):: from_fields
-  !! Resulting extrapolated field. May be the same field or a field on
-  !! a different mesh (different degree).
-  type(scalar_field), dimension(:), intent(inout):: to_fields
-  !! positions, and upward normal vector on the whole domain
-  type(vector_field), target, intent(inout):: positions
-  type(vector_field), target, intent(in):: vertical_normal
-
-  !! the surface elements (faces numbers) that make up the surface
-  integer, dimension(:), intent(in):: surface_element_list
-  !! If provided the projected surface mesh onto horizontal coordinates 
-  !! and its associated rtree/pickers are cached under this name and 
-  !! attached to the 'positions'. In this case when called again with
-  !! the same 'positions' and the same surface_name,
-  !! the same surface_element_list should again be provided.
-  character(len=*), optional, intent(in):: surface_name
-
-  character(len=FIELD_NAME_LEN):: lsurface_name
-  real, dimension(:,:), allocatable:: loc_coords
-  integer, dimension(:), allocatable:: seles
-  integer i, j, to_nodes, face
-
-  assert(size(from_fields)==size(to_fields))
-  assert(element_count(from_fields(1))==element_count(to_fields(1)))
-  do i=2, size(to_fields)
-     assert(to_fields(1)%mesh==to_fields(i)%mesh)
-     assert(from_fields(1)%mesh==from_fields(i)%mesh)
-  end do
-  
-  to_nodes=nowned_nodes(to_fields(1))
-  ! local coordinates is one more than horizontal coordinate dim
-  allocate( seles(to_nodes), loc_coords(1:positions%dim, 1:to_nodes) )
-  
-  if (present(surface_name)) then
-    lsurface_name=surface_name
-  else
-    lsurface_name="TempSurfaceName"
-  end if
-  
-  ! project the positions of to_fields(1)%mesh into the horizontal plane
-  ! and returns 'seles' (indices in surface_element_list) and loc_coords 
-  ! to tell where these projected nodes are found in the surface mesh
-  call horizontal_picker(to_fields(1)%mesh, positions, vertical_normal, &
-      surface_element_list, lsurface_name, &
-      seles, loc_coords)
-  
-  ! interpolate using the returned faces and loc_coords
-  do i=1, size(to_fields)
-    do j=1, to_nodes
-      face=surface_element_list(seles(j))
-      call set(to_fields(i), j, &
-           dot_product( eval_shape( face_shape(from_fields(i), face), loc_coords(:,j) ), &
-             face_val( from_fields(i), face ) ))
-     end do
-  end do
-    
-  if (IsParallel()) then
-    do i=1, size(to_fields)
-      call halo_update(to_fields(i))
-    end do
-  end if
-  
-  if (.not. present(surface_name)) then
-    call remove_boundary_condition(positions, "TempSurfaceName")
-  end if
-  
-end subroutine VerticalExtrapolationMultiple
-  
-subroutine horizontal_picker(mesh, positions, vertical_normal, &
-  surface_element_list, surface_name, &
-  seles, loc_coords)
-  !! Searches the nodes of 'mesh' in the surface mesh above.
-  !! Returns the surface elements 'seles' that each node lies under
-  !! and the loc_coords in this element of this node projected
-  !! upward (radially on the sphere) onto the surface mesh
-  
-  !! mesh
-  type(mesh_type), intent(in):: mesh
-  !! a valid positions field for the whole domain, not necessarily on 'mesh'
-  !! for instance in a periodic domain, mesh is periodic and positions should not be
-  type(vector_field), target, intent(inout):: positions
-  !! upward normal vector on the whole domain
-  type(vector_field), target, intent(in):: vertical_normal
-  !! the surface elements (faces numbers) that make up the surface
-  integer, dimension(:), intent(in):: surface_element_list
-  !! The projected surface mesh onto horizontal coordinates 
-  !! and its associated rtree/pickers are cached under this name and 
-  !! attached to the 'positions'. When called again with
-  !! the same 'positions' and the same surface_name,
-  !! the same surface_element_list should again be provided.
-  character(len=*), intent(in):: surface_name
-  
-  !! returned surface elements (face numbers in 'mesh')
-  !! and loc coords each node has been found in
-  !! size(seles)==size(loc_coords,2)==nowned_nodes(mesh)
-  integer, dimension(:), intent(out):: seles
-  real, dimension(:,:), intent(out):: loc_coords
-  
-  type(vector_field):: mesh_positions
-  type(vector_field), pointer:: horizontal_positions
-  integer, dimension(:), pointer:: horizontal_mesh_list
-  real, dimension(:,:), allocatable:: horizontal_coordinate
-  real, dimension(vertical_normal%dim):: normal_vector
-  real, dimension(positions%dim):: xyz
-  integer:: i, stat, nodes
-  
-  assert(.not. mesh_periodic(positions))
-  
-  ! search only the first owned nodes
-  nodes=nowned_nodes(mesh)
-  
-  assert( size(seles)==nodes )
-  assert(size(loc_coords,1)==positions%dim)
-  assert( size(loc_coords,2)==nodes )
-  
-  if (mesh==positions%mesh) then
-    mesh_positions=positions
-    ! make mesh_positions indep. ref. of the field, so we can deallocate it 
-    ! safely without destroying positions
-    call incref(mesh_positions)
-  else
-    call allocate(mesh_positions, positions%dim, mesh, &
-      name='ToPositions_VerticalExtrapolation')
-    call remap_field(positions, mesh_positions, stat)
-    if (stat/=0 .and. stat/=REMAP_ERR_HIGHER_LOWER_CONTINUOUS .and. &
-      stat/=REMAP_ERR_UNPERIODIC_PERIODIC) then
-      ! Mapping from higher order to lower order is allowed for coordinates
-      ! (well depends on how the higher order is derived from the lower order)
-      !
-      ! Mapping to periodic coordinates is ok in this case as we only need
-      ! locations for the nodes individually (i.e. we don't care about elements
-      ! in 'mesh') - the created horizontal_positions will be non-periodic
-      ! So that we should be able to find nodes on the periodic boundary on
-      ! either side. Using this to interpolate is consistent as long as
-      ! the interpolated from field is indeed periodic.
-      FLAbort("Unknown error in remmaping positions in horizontal_picker.")
-    end if
-  end if
-  
-  
-  ! create an array of the horizontally projected coordinates of the 'mesh' nodes
-  allocate( horizontal_coordinate(1:positions%dim-1, 1:nodes) )
-  if (have_option('/geometry/spherical_earth')) then
-    assert(mesh_positions%dim==3)
-    do i=1, nodes
-      xyz=node_val(mesh_positions, i)
-      if (xyz(3)>0.0) then
-        horizontal_coordinate(:,i) = &
-          map2horizontal_sphere( node_val(mesh_positions, i), +1.0)
-      else
-        horizontal_coordinate(:,i) = &
-          map2horizontal_sphere( node_val(mesh_positions, i), -1.0)
-        horizontal_coordinate(1,i)=horizontal_coordinate(1,i)+HEMISPHERE_SHIFT
-      end if
-    end do
-  else
-    assert( vertical_normal%field_type==FIELD_TYPE_CONSTANT )
-    normal_vector=node_val(vertical_normal,1)
-
-    do i=1, nodes
-      horizontal_coordinate(:,i) = &
-          map2horizontal( node_val(mesh_positions, i), normal_vector )
-    end do
-  end if
-    
-  call get_horizontal_positions(positions, surface_element_list, &
-      vertical_normal, surface_name, &
-      horizontal_positions, horizontal_mesh_list)
-    
-  call picker_inquire( horizontal_positions, horizontal_coordinate, &
-    seles, loc_coords, global=.false. )
-    
-  ! in the spherical case some of the surface elements may be duplicated
-  ! within horizontal positions, the returned seles should refer to entries
-  ! in surface_element_list however - also check for nodes not found
-  do i=1, size(seles)
-    if (seles(i)>0) then
-      seles(i)=horizontal_mesh_list(seles(i))
-    else
-      ewrite(-1,*) "For node with coordinate", node_val(mesh_positions, i)
-      ewrite(-1,*) "no top surface node was found."
-      FLAbort("Something wrong with the geometry.")
-    end if
-  end do
-    
-  call deallocate(mesh_positions)
-  deallocate(horizontal_mesh_list)
-
-end subroutine horizontal_picker
-  
-subroutine get_horizontal_positions(positions, surface_element_list, vertical_normal, surface_name, &
-  horizontal_positions, horizontal_mesh_list)
-! returns a horizontal positions field over the surface mesh indicated by
-! 'surface_element_list'. This field will be created and cached on 'positions'
-! as a surface field attached to a dummy boundary condition under the name surface_name
-  type(vector_field), intent(inout):: positions
-  type(vector_field), intent(in):: vertical_normal
-  integer, dimension(:), intent(in):: surface_element_list
-  character(len=*), intent(in):: surface_name
-  
-  ! Returns the horizontal positions on a horizontal mesh, this mesh
-  ! may have some of the faces in surface_element_list duplicated -
-  ! this is only in the spherical case where the horizontal mesh
-  ! consists of two disjoint regions representing the two hemispheres
-  ! and surface elements near the equator may be represented in both.
-  ! Therefore we also return a map between elements in the horizontal positions mesh
-  ! and entries in surface_element_list. If ele is an element in the horizontal 
-  ! position mesh, then surface_element_list(horizontal_mesh(ele)) 
-  ! is the face number in the full  mesh
-  ! horizontal_positions is a borrowed reference, don't allocate
-  ! horizontal_mesh_list does need to be deallocated
-  type(vector_field), pointer:: horizontal_positions
-  integer, dimension(:), pointer:: horizontal_mesh_list
-  
-  integer, dimension(:), pointer:: horizontal_sele_list
-  integer:: i, j, k, sele
-  
-  if (.not. has_boundary_condition_name(positions, surface_name)) then
-    if (have_option('/geometry/spherical_earth')) then
-      call create_horizontal_positions_sphere(positions, &
-        surface_element_list, surface_name)
-    else
-      call create_horizontal_positions_flat(positions, &
-        surface_element_list, vertical_normal, surface_name)
-    end if
-  end if
-    
-  horizontal_positions => extract_surface_field(positions, surface_name, &
-    trim(surface_name)//"HorizontalCoordinate")
-    
-!   call vtk_write_fields('horizontal_mesh', 0,      horizontal_positions, &
-!     horizontal_positions%mesh)    
-    
-    
-  allocate( horizontal_mesh_list(1:element_count(horizontal_positions)) )
-  if (have_option('/geometry/spherical_earth')) then
-    call get_boundary_condition(positions, surface_name, &
-      surface_element_list=horizontal_sele_list)
-    ! by construction the map can be obtained by running 
-    ! through surface_element_list twice (for each hemisphere)
-    i=1 ! position in horizontal_sele_list
-    sele=horizontal_sele_list(i) ! face we're looking for
-outer_loop: &
-    do j=1, 2
-      do k=1, size(surface_element_list)
-        if (surface_element_list(k)==sele) then
-          ! found matching position in surface_element_list
-          horizontal_mesh_list(i)=k
-          ! next one to search
-          i=i+1
-          if (i>size(horizontal_sele_list)) exit outer_loop
-          sele=horizontal_sele_list(i)
-        end if
-      end do
-    end do outer_loop
-      
-    if (i<=size(horizontal_sele_list)) then
-      ! not all were found in 2 loops through surface_element_list
-      ! something's wrong
-      FLAbort("Internal error in horizontal mesh administration")
-    end if
-        
-  else
-    ! no duplication: horizontal_mesh_list is simply the identity map
-    do i=1, size(horizontal_mesh_list)
-      horizontal_mesh_list(i)=i
-    end do
-  end if
-  
-end subroutine get_horizontal_positions
-  
-subroutine create_horizontal_positions_flat(positions, surface_element_list, vertical_normal, surface_name)
-! adds a "boundary condition" to 'positions' with an associated vector surface field containing a dim-1 
-! horizontal coordinate field that can be used to map from the surface mesh specified 
-! by 'surface_element_list'. This "boundary condition" will be stored under the name 'surface_name'
-! The horizontal coordinates are created by projecting out the component
-! in the direction of 'vertical_normal' and then throwing out the x, y or z 
-! coordinate that is most aligned with 'vertical_normal'
-  type(vector_field), intent(inout):: positions
-  type(vector_field), intent(in):: vertical_normal
-  integer, dimension(:), intent(in):: surface_element_list
-  character(len=*), intent(in):: surface_name
-    
-  type(mesh_type), pointer:: surface_mesh
-  type(vector_field):: horizontal_positions
-  integer, dimension(:), pointer:: surface_node_list
-  real, dimension(vertical_normal%dim):: normal_vector
-  integer:: i, node
-  
-  assert(vertical_normal%field_type==FIELD_TYPE_CONSTANT)
-  normal_vector=node_val(vertical_normal, 1)
-  
-  call add_boundary_condition_surface_elements(positions, &
-    name=surface_name, type="verticalextrapolation", &
-    surface_element_list=surface_element_list)
-  ! now get back the created surface mesh
-  ! and surface_node_list a mapping between node nos in the projected mesh and node nos in the original 'positions' mesh
-  call get_boundary_condition(positions, name=surface_name, &
-    surface_mesh=surface_mesh, surface_node_list=surface_node_list)
-  call allocate( horizontal_positions, dim=positions%dim-1, mesh=surface_mesh, &
-    name=trim(surface_name)//"HorizontalCoordinate" )
-    
-  call insert_surface_field(positions, name=surface_name, &
-    surface_field=horizontal_positions)
-    
-  do i=1, size(surface_node_list)
-    node=surface_node_list(i)
-    call set(horizontal_positions, i, &
-        map2horizontal(node_val(positions, node), normal_vector))
-  end do
-    
-  call deallocate( horizontal_positions )
-
-end subroutine create_horizontal_positions_flat
-
-subroutine create_horizontal_positions_sphere(positions, surface_element_list, surface_name)
-! adds a "boundary condition" to 'positions' with an associated vector surface field containing a dim-1 
-! horizontal coordinate field that can be used to map from the surface mesh specified 
-! by 'surface_element_list'. This "boundary condition" will be stored under the name 'surface_name'
-! The horizontal coordinates are created by a stereographic projection from the unit sphere to the
-! xy-plane. The northern hemisphere is projected to the unit-circle (including some extra surface
-! elements on the southern hemisphere around the equator). The southern hemisphere (including some
-! northern surface elements near the equator) is also projected to a unit circle but translated 
-! away from the origin to not overlap with the projected northern hemisphere. Thus the created 
-! projected horizontal positions field consists of two disjoint areas in the plane corresponding 
-! to both hemispheres, and elements in the original surface mesh (near the equator) may appear 
-! twice in the projected field.
-  type(vector_field), intent(inout):: positions
-  integer, dimension(:), intent(in):: surface_element_list
-  character(len=*), intent(in):: surface_name
-  
-  type(vector_field), pointer:: horizontal_positions_north, horizontal_positions_south
-  type(vector_field):: horizontal_positions
-  type(mesh_type), pointer:: surface_mesh_north, surface_mesh_south
-  type(mesh_type):: surface_mesh
-  integer, dimension(:), pointer:: surface_element_list_north, surface_element_list_south
-  integer, dimension(:), allocatable:: surface_element_list_combined
-  integer:: nodes_north, elements_south, elements_north
-  integer:: i
-  
-  ! first create 2 separate horizontal coordinate fields
-  ! for each hemisphere
-  
-  call create_horizontal_positions_hemisphere(positions, &
-    surface_element_list, trim(surface_name)//'North', +1.0)
-    
-  call create_horizontal_positions_hemisphere(positions, &
-    surface_element_list, trim(surface_name)//'South', -1.0)
-  
-  ! these are stored as bcs under the positions
-  ! retrieve this information  back:
-  
-  call get_boundary_condition(positions, name=trim(surface_name)//'North', &
-    surface_mesh=surface_mesh_north, &
-    surface_element_list=surface_element_list_north)
-  call get_boundary_condition(positions, name=trim(surface_name)//'South', &
-    surface_mesh=surface_mesh_south, &
-    surface_element_list=surface_element_list_south)    
-    
-  horizontal_positions_north => extract_surface_field(positions, &
-      trim(surface_name)//'North', trim(surface_name)//"NorthHorizontalCoordinate")
-  horizontal_positions_south => extract_surface_field(positions, &
-      trim(surface_name)//'South', trim(surface_name)//"SouthHorizontalCoordinate")    
-    
-  ! merge these 2 meshes
-  surface_mesh=merge_meshes( (/ surface_mesh_north, surface_mesh_south /) )
-  
-  ! and merge the positions field
-  call allocate( horizontal_positions, positions%dim-1, surface_mesh, &
-    trim(surface_name)//"HorizontalCoordinate" )
-    
-  nodes_north=node_count(surface_mesh_north)
-  do i=1, nodes_north
-    call set( horizontal_positions, i, &
-      node_val(horizontal_positions_north, i))
-  end do
-  ! but translate the projected southern hemisphere to the left
-  ! to not overlap it with the northern hemisphere
-  do i=1, node_count(surface_mesh_south)
-    call set( horizontal_positions, nodes_north+i, &
-      node_val(horizontal_positions_south, i)+(/ HEMISPHERE_SHIFT, 0.0 /) )
-  end do
-    
-  ! merge the surface element lists
-  ! (note that this is different than the incoming surface_element_list
-  ! as it will have some duplicate equatorial elements)
-  elements_north=size(surface_element_list_north)
-  elements_south=size(surface_element_list_south)
-  allocate(surface_element_list_combined(1:elements_north+elements_south))
-  surface_element_list_combined(1:elements_north)=surface_element_list_north
-  surface_element_list_combined(elements_north+1:)=surface_element_list_south
-  
-  ! finally the bc for the combined surface mesh:
-  ! (note that this creates a new surface mesh different than
-  ! the merged mesh, that we won't be using)
-  call add_boundary_condition_surface_elements( &
-    positions, name=surface_name, type="verticalextrapolation", &
-    surface_element_list=surface_element_list_combined)
-  
-  ! insert the horizontal positions under this bc
-  call insert_surface_field( positions, name=surface_name, &
-    surface_field=horizontal_positions)
-  
-  ! everything is safely stored, so we can deallocate our references
-  call deallocate(horizontal_positions)
-  call deallocate(surface_mesh)
-  deallocate(surface_element_list_combined)
-  
-  ! also we won't need the 2 hemisphere bcs anymore
-  call remove_boundary_condition( positions, trim(surface_name)//'North')
-  call remove_boundary_condition( positions, trim(surface_name)//'South')
-    
-end subroutine create_horizontal_positions_sphere
-
-subroutine create_horizontal_positions_hemisphere(positions, &
-  surface_element_list, surface_name, hemi_sign)
-  type(vector_field), intent(inout):: positions
-  integer, dimension(:), intent(in):: surface_element_list
-  character(len=*), intent(in):: surface_name
-  real, intent(in):: hemi_sign
-
-  type(vector_field):: horizontal_positions
-  type(mesh_type), pointer:: surface_mesh
-  real, dimension(2):: xy
-  integer, dimension(:), pointer:: nodes, surface_node_list
-  integer:: i, j, sele, node
-  
-  type(integer_set):: surface_element_set
-  
-  call allocate(surface_element_set)
-  
-  do i=1, size(surface_element_list)
-    sele=surface_element_list(i)
-    if (sele_is_on_hemisphere(sele)) then
-      call insert(surface_element_set, sele)
-    end if
-  end do
-    
-  call add_boundary_condition_surface_elements(positions, &
-    name=surface_name, type="verticalextrapolation", &
-    surface_element_list=set2vector(surface_element_set))
-  call deallocate(surface_element_set)
-  
-  call get_boundary_condition(positions, name=surface_name, &
-    surface_mesh=surface_mesh, surface_node_list=surface_node_list)
-    
-  call allocate( horizontal_positions, dim=2, mesh=surface_mesh, &
-    name=trim(surface_name)//"HorizontalCoordinate" )
-    
-  call insert_surface_field(positions, name=surface_name, &
-    surface_field=horizontal_positions)
-    
-  do i=1, element_count(horizontal_positions)
-    nodes => ele_nodes(horizontal_positions, i)
-    do j=1, size(nodes)
-      node=surface_node_list( nodes(j) )
-      xy=map2horizontal_sphere(node_val(positions, node), hemi_sign)
-      assert(abs(xy(1))<HEMISPHERE_SHIFT/2.0)
-      call set( horizontal_positions, nodes(j), xy)
-    end do
-  end do
-    
-  call deallocate(horizontal_positions)
-
-  contains
-
-  logical function sele_is_on_hemisphere(sele)
-    integer, intent(in):: sele
-
-    real, dimension(face_loc(positions, sele)):: znodes
-
-    znodes=face_val(positions, 3, sele)
-    ! if any of the nodes is on the hemisphere or *very* close to the equator
-    ! (relative to the other nodes)    
-    sele_is_on_hemisphere = any(hemi_sign * znodes > -sum(abs(znodes))*1e-6)
-
-  end function sele_is_on_hemisphere
-
-end subroutine create_horizontal_positions_hemisphere
-
-function map2horizontal(xyz, normal_vector)
-real, dimension(:), intent(in):: xyz, normal_vector
-real, dimension(size(xyz)-1):: map2horizontal
-
-  real, dimension(size(xyz)):: hxyz
-  integer:: i, c, takeout
-  
-  ! first subtract of the vertical component
-  hxyz=xyz-dot_product(xyz, normal_vector)*normal_vector
-  
-  ! then leave out the "most vertical" coordinate
-  takeout=maxloc(abs(normal_vector), dim=1)
-  
-  c=1
-  do i=1, size(xyz)
-    if (i==takeout) cycle
-    map2horizontal(c)=hxyz(i)
-    c=c+1
-  end do
-
-end function map2horizontal
-
-function map2horizontal_sphere(xyz, hemi_sign)
-real, dimension(3), intent(in):: xyz
-real, intent(in):: hemi_sign
-real, dimension(2):: map2horizontal_sphere
-
-  real:: r
-  
-  r=sqrt(sum(xyz**2))
-  map2horizontal_sphere=xyz(1:2)/(r+hemi_sign*xyz(3))
-    
-end function map2horizontal_sphere
-
-function VerticalProlongationOperator(mesh, positions, vertical_normal, &
-  surface_element_list, surface_mesh)
-  !! creates a prolongation operator that prolongates values on
-  !! a surface mesh to a full mesh below using the same interpolation
-  !! as the vertical extrapolation code above. The transpose of this prolongation
-  !! operator can be used as a restriction/clustering operator 
-  !! from the full mesh to the surface.
-  type(csr_matrix) :: VerticalProlongationOperator
-  !! the mesh to which to prolongate, its nodes are only considered as
-  !! a set of loose points
-  type(mesh_type), target, intent(in):: mesh
-  !! positions on the whole domain (doesn't have to be the same mesh)
-  type(vector_field), intent(inout):: positions
-  !! upward normal vector on the whole domain
-  type(vector_field), target, intent(in):: vertical_normal
-  !! the face numbers of the surface mesh
-  integer, dimension(:), intent(in):: surface_element_list
-  !! Optionally a surface mesh may be provided that represents the nodes
-  !! /from/ which to interpolate (may be of different signature than 'mesh')
-  !! Each column in the prolongator will correspond to a node in this surface_mesh.
-  !! If not provided, the "from mesh" is considered to consist of faces
-  !! given by surface_element_list (and same cont. and order as 'mesh').
-  !! In this case however empty columns, i.e. surface nodes from which no
-  !! value in the full mesh is interpolated, are removed and there is no
-  !! necessary relation between column numbering and surface node numbering.
-  type(mesh_type), intent(in), target, optional:: surface_mesh
-
-  type(csr_sparsity):: sparsity
-  type(mesh_type), pointer:: lsurface_mesh
-  real, dimension(:,:), allocatable:: loc_coords
-  real, dimension(:), allocatable:: mat
-  real:: coef
-  integer, dimension(:), pointer:: snodes
-  integer, dimension(:), allocatable:: seles, colm, findrm, snod2used_snod
-  integer i, j, k, rows, entries, count, snod, sele
-  ! coefficient have to be at least this otherwise they're negligable in an interpolation
-  real, parameter :: COEF_EPS=1d-10
-  
-  ! only assemble the rows associted with nodes we own
-  rows=nowned_nodes(mesh)
-  
-  ! local coordinates is one more than horizontal coordinate dim
-  allocate( seles(rows), loc_coords(1:positions%dim, 1:rows) )
-  
-  ! project the positions of to_fields(1)%mesh into the horizontal plane
-  ! and returns 'seles' (indices in surface_element_list) and loc_coords 
-  ! to tell where these projected nodes are found in the surface mesh
-  call horizontal_picker(mesh, positions, vertical_normal, &
-      surface_element_list, "TempSurfaceName", &
-      seles, loc_coords)
-
-  ! count upper estimate for n/o entries for sparsity
-  entries=0
-  do i=1, rows
-     entries=entries+face_loc(mesh, seles(i))
-  end do
-  ! preliminary matrix:  
-  allocate( mat(1:entries), findrm(1:rows+1), &
-    colm(1:entries) )
-  
-  if (.not. present(surface_mesh)) then
-     ! We use the entire surface mesh of 'mesh'
-     lsurface_mesh => mesh%faces%surface_mesh
-     ! Not all surface nodes may be used (i.e. interpolated from) - even
-     ! within surface elements that /are/ in surface_element-list.
-     ! We need a map between global surface node numbering and 
-     ! a consecutive numbering of used surface nodes.
-     ! (this will be the column numbering)
-     allocate(snod2used_snod(1:node_count(lsurface_mesh)))
-     snod2used_snod=0
-     count=0 ! counts number of used surface nodes
-  else
-     lsurface_mesh => surface_mesh
-  end if
-
-  entries=0 ! this time only count nonzero entries
-  do i=1, rows
-
-     ! beginning of each row in mat
-     findrm(i)=entries+1
-      
-     if (present(surface_mesh)) then
-       ! element number within surface_mesh
-       sele=seles(i)
-     else
-       ! face number in 'mesh', i.e. element number within entire surface_mesh
-       sele=surface_element_list(seles(i))
-     end if
-     
-     snodes => ele_nodes(lsurface_mesh, sele)
-     
-     do j=1, size(snodes)
-       coef=eval_shape(ele_shape(lsurface_mesh, sele), j, loc_coords(:,i))
-       snod=snodes(j)
-       if (abs(coef)>COEF_EPS) then
-         if (.not. present(surface_mesh)) then
-           if (snod2used_snod(snod)==0) then
-             ! as of yet unused surface node
-             count=count+1
-             snod2used_snod(snod)=count
-           end if
-           ! this is the column index we're gonna use instead
-           snod=snod2used_snod(snod)
-         end if
-         entries=entries+1
-         colm(entries)=snod
-         mat(entries)=coef
-       end if
-     end do
-       
-  end do
-  findrm(i)=entries+1
-  
-  if (present(surface_mesh)) then
-    ! we haven't counted used surface nodes, instead we're using all 
-    ! nodes of surface mesh as columns
-    count=node_count(surface_mesh)
-  end if
-  
-  call allocate(sparsity, rows, count, &
-     entries, diag=.false., name="VerticalProlongationSparsity")
-  sparsity%findrm=findrm
-  sparsity%colm=colm(1:entries)
-  ! for lots of applications it's good to have sorted rows
-  call sparsity_sort(sparsity)
-  
-  call allocate(VerticalProlongationOperator, sparsity, &
-    name="VerticalProlongationOperator")
-  call deallocate(sparsity)
-
-  ! as the sparsity has been sorted the ordering of mat(:) no longer
-  ! matches that of sparsity%colm, however it still matches the original 
-  ! unsorted colm(:)
-  do i=1, rows
-     do k=findrm(i), findrm(i+1)-1
-       j=colm(k)
-       call set(VerticalProlongationOperator, i, j, mat(k))
-     end do
-  end do
-  
-  if (.not. present(surface_mesh)) then
-    deallocate( snod2used_snod )
-  end if
-  deallocate( findrm, colm, mat )
-  deallocate( seles, loc_coords )
-  
-  call remove_boundary_condition(positions, "TempSurfaceName")
-  
-end function VerticalProlongationOperator
-  
-subroutine vertical_element_ordering(ordered_elements, face_normal_gravity, optimal_ordering)
-!!< Calculates an element ordering such that each element is
-!!< is preceded by all elements above it.
-integer, dimension(:), intent(out):: ordered_elements
-!! need to supply face_normal_gravity matrix,
-!! created by compute_face_normal_gravity() subroutine below
-type(csr_matrix), intent(in):: face_normal_gravity
-!! returns .true. if an optimal ordering is found, i.e there are no 
-!! cycles, i.o.w. elements that are (indirectly) above and below each other 
-!! at the same time (deck of cards problem).
-logical, optional, intent(out):: optimal_ordering
-  
-  type(dynamic_bin_type) dbin
-  real, dimension(:), pointer:: inn
-  integer, dimension(:), pointer:: neigh
-  integer, dimension(:), allocatable:: bin_list
-  integer i, j, elm, bin_no
-  logical warning
-  
-  assert( size(ordered_elements)==size(face_normal_gravity,1) )
-  
-  ! create binlist, i.e. assign each element to a bin, according to
-  ! the number of elements above it
-  allocate(bin_list(1:size(ordered_elements)))
-  do i=1, size(ordered_elements)
-    neigh => row_m_ptr(face_normal_gravity, i)
-    inn => row_val_ptr(face_normal_gravity, i)
-    ! elements with no element above it go in bin 1
-    ! elements with n elements above it go in bin n+1
-    ! neigh>0 so we don't count exterior boundary faces
-    bin_list(i)=count( inn<-VERTICAL_INTEGRATION_EPS .and. neigh>0 )+1
-  end do
-    
-  call allocate(dbin, bin_list)
-  
-  warning=.false.
-  do i=1, size(ordered_elements)
-    ! pull an element from the first non-empty bin
-    ! (hopefully an element with no unprocessed elements above)
-    call pull_element(dbin, elm, bin_no)
-    ordered_elements(i)=elm
-    ! if this is bin one then it is indeed an element with no unprocessed 
-    !  elements above, otherwise issue a warning
-    if (bin_no>1) warning=.true.
-    
-    ! update elements below:
-    
-    ! adjacent elements:
-    neigh => row_m_ptr(face_normal_gravity, elm)
-    inn => row_val_ptr(face_normal_gravity, elm)
-    do j=1, size(neigh)
-       if (inn(j)>VERTICAL_INTEGRATION_EPS .and. neigh(j)>0) then
-         ! element neigh(j) is below element i, therefore now has one
-         ! less unprocessed element above it, so can be moved to
-         ! lower bin.
-         if (.not. element_pulled(dbin, neigh(j))) then
-           ! but only if neigh(j) itself hasn't been selected yet
-           ! (which might happen for imperfect vertical orderings)
-           call move_element(dbin, neigh(j), bin_list(neigh(j))-1)
-         end if
-       end if
-    end do
-  end do
-  
-  if (warning) then
-    ! this warning may be reduced (in verbosity level) if it occurs frequently:
-    ewrite(-1,*) "Warning: vertical_element_ordering has detected a cycle."
-    ewrite(-1,*) "(deck of cards problem). This may reduce the efficiency"
-    ewrite(-1,*) "of your vertically sweeping solve."
-  end if
-  
-  if (present(optimal_ordering)) then
-    optimal_ordering=.not. warning
-  end if
-  
-  call deallocate(dbin)
-  
-end subroutine vertical_element_ordering
-  
-subroutine compute_face_normal_gravity(face_normal_gravity, &
-  positions, vertical_normal)
-!!< Returns a matrix where A_ij is the inner product of the face normal
-!!< and the gravity normal vector of the face between element i and j.
-type(csr_matrix), intent(out):: face_normal_gravity
-type(vector_field), target, intent(in):: positions, vertical_normal
-
-  type(mesh_type), pointer:: mesh
-  real, dimension(:), pointer:: face_normal_gravity_val
-  real, dimension(:), allocatable:: detwei_f
-  real, dimension(:,:), allocatable:: face_normal, gravity_normal
-  integer, dimension(:), pointer:: neigh, faces
-  real inn, area
-  integer sngi, nloc, i, k
-    
-  mesh => positions%mesh
-  call allocate(face_normal_gravity, mesh%faces%face_list%sparsity)
-  call zero(face_normal_gravity)
-  
-  sngi=face_ngi(mesh, 1)
-  nloc=ele_loc(mesh,1)
-  allocate( detwei_f(1:sngi), &
-    face_normal(1:positions%dim, 1:sngi), &
-    gravity_normal(1:positions%dim, 1:sngi))
-  
-  do i=1, element_count(mesh)
-     ! elements adjacent to element i
-     ! this is a row (column indices) in the mesh%faces%face_list matrix
-     neigh => ele_neigh(mesh, i)
-     ! the surrounding faces
-     ! this is a row (integer values) in the mesh%faces%face_list matrix
-     faces => ele_faces(mesh, i)
-     do k=1, size(neigh)
-        if (neigh(k)>i .or. neigh(k)<=0) then
-           ! only handling neigh(k)>i to ensure anti-symmetry of the matrix
-           ! (and more efficient of course)
-           call transform_facet_to_physical(positions, faces(k), &
-                detwei_f=detwei_f, &
-                normal=face_normal)
-           gravity_normal=face_val_at_quad(vertical_normal, faces(k))
-           area=sum(detwei_f)
-           ! inner product of face normal and vertical normal
-           ! integrated over face
-           inn=sum(matmul(face_normal*gravity_normal, detwei_f))/area
-           if (neigh(k)>0) then
-              call set(face_normal_gravity, i, neigh(k), inn)
-              call set(face_normal_gravity, neigh(k), i, -inn)
-           else
-              ! exterior surface: matrix entry does not have valid
-              ! column index, still want to store its value, so we
-              ! use a pointer
-              face_normal_gravity_val => row_val_ptr(face_normal_gravity, i)
-              face_normal_gravity_val(k)=inn
-           end if
-        end if
-     end do
-        
-  end do
-
-end subroutine compute_face_normal_gravity
-  
-subroutine vertical_integration_scalar(from_field, to_field, &
-    positions, vertical_normal, surface_element_list, rhs)
-!!< See description vertical_integration_multiple
-type(scalar_field), intent(in):: from_field
-type(scalar_field), intent(in):: to_field
-type(vector_field), intent(in):: positions, vertical_normal
-integer, dimension(:), intent(in):: surface_element_list
-type(scalar_field), optional, intent(in):: rhs
-
-  type(scalar_field) to_fields(1)
-    
-  to_fields=(/ to_field /)
-  if (present(rhs)) then     
-     call vertical_integration_multiple( (/ from_field /), to_fields, &
-        positions, vertical_normal, surface_element_list, rhs=(/ rhs /) )
-  else  
-     call vertical_integration_multiple( (/ from_field /), to_fields, &
-        positions, vertical_normal, surface_element_list)
-  end if
-
-end subroutine vertical_integration_scalar
-
-subroutine vertical_integration_vector(from_field, to_field, &
-    positions, vertical_normal, surface_element_list, rhs)
-!!< See description vertical_integration_multiple
-type(vector_field), intent(in):: from_field
-type(vector_field), intent(in):: to_field
-type(vector_field), intent(in):: positions, vertical_normal
-integer, dimension(:), intent(in):: surface_element_list
-type(vector_field), optional, intent(in):: rhs
-  
-  type(scalar_field), dimension(from_field%dim):: from_field_components, &
-     to_field_components, rhs_components
-  integer i
-  
-  assert(from_field%dim==to_field%dim)
-  
-  do i=1, from_field%dim
-     from_field_components(i)=extract_scalar_field(from_field, i)
-     to_field_components(i)=extract_scalar_field(to_field, i)
-     if (present(rhs)) then
-        rhs_components(i)=extract_scalar_field(rhs, i)
-     end if
-  end do
-  
-  if (present(rhs)) then
-     call vertical_integration_multiple( from_field_components, &
-        to_field_components, positions, vertical_normal, &
-        surface_element_list, rhs=rhs_components)
-  else
-     call vertical_integration_multiple( from_field_components, &
-        to_field_components, positions, vertical_normal, &
-        surface_element_list)
-  end if
-  
-end subroutine vertical_integration_vector
-
-subroutine vertical_integration_multiple(from_fields, to_fields, &
-    positions, vertical_normal, surface_element_list, rhs)
-!!< This subroutine solves: dP/dz=rhs using DG
-!!< It can be used for vertical integration downwards (dP/dz=0) as a drop
-!!< in replacement of VerticalExtrapolation hence its similar interface.
-!!< The field P is the to_field. A boundary condition is given by the 
-!!< from_field. Again completely similar to VerticalExtrapolation, it may
-!!< be defined as a surface field on surface elements given by 
-!!< surface_element_list or it may be a field on the complete mesh
-!!< in which case only its values on these surface elements are used.
-!!< vertical_normal specifies the direction in which to integrate (usually downwards)
-!!< If not specified rhs is assumed zero.
-!!<
-!!< This version accepts multiple from_fields, to_fields and rhs
-type(scalar_field), dimension(:), intent(in):: from_fields
-type(scalar_field), dimension(:), intent(inout):: to_fields
-type(vector_field), intent(in):: positions, vertical_normal
-integer, dimension(:), intent(in):: surface_element_list
-type(scalar_field), dimension(:), optional, intent(in):: rhs
-  
-  type(csr_matrix) face_normal_gravity
-  type(element_type), pointer:: ele_shp, face_shp, x_face_shp
-  real, dimension(:), pointer:: inn
-  real, dimension(:,:,:), allocatable:: surface_rhs, dele_shp
-  real, dimension(:,:), allocatable:: ele_mat, face_mat, ele_rhs
-  real, dimension(:), allocatable:: detwei, detwei_f
-  integer, dimension(:), pointer:: neigh, ele_nds, faces, face_lnds
-  integer, dimension(:), allocatable:: ordered_elements, face_nds, face_nds2
-  integer nloc, snloc, ngi, sngi
-  integer i, j, k, f, f2, elm, it, noit
-  logical optimal_ordering, from_surface_fields
-  
-  assert( size(from_fields)==size(to_fields) )
-  
-  ! computes inner product of face normal and gravity (see above)
-  call compute_face_normal_gravity(face_normal_gravity, &
-     positions, vertical_normal)
-  
-  ! determine an ordering for the elements based on this
-  allocate( ordered_elements(1:element_count(positions)) )
-  call vertical_element_ordering(ordered_elements, face_normal_gravity, &
-     optimal_ordering)
-     
-  ! General initalisation
-  !-----------------------
-  ! various grid numbers
-  nloc=ele_loc(to_fields(1), 1)
-  snloc=face_loc(to_fields(1), 1)
-  ngi=ele_ngi(positions, 1)
-  sngi=face_ngi(positions, 1)
-  ! shape functions
-  ele_shp => ele_shape(to_fields(1), 1)
-  face_shp => face_shape(to_fields(1), 1)
-  x_face_shp => face_shape(positions, 1)
-  ! various allocations:
-  allocate( &
-     surface_rhs(1:snloc, 1:size(to_fields), 1:surface_element_count(positions)), &
-     ele_mat(1:nloc, 1:nloc), ele_rhs(1:nloc,1:size(to_fields)), &
-     dele_shp(1:nloc, 1:ngi, 1:positions%dim), &
-     face_mat(1:snloc, 1:snloc), face_nds(1:snloc), face_nds2(1:snloc), &
-     detwei(1:ngi), detwei_f(1:sngi))
-     
-  if (element_count(from_fields(1))==size(surface_element_list)) then
-     ! from_fields are fields over the surface mesh only
-     ! so we're using all of its values:
-     from_surface_fields=.true.
-     ! check the other fields as well:
-     do k=2, size(from_fields)
-       assert( element_count(from_fields(k))==size(surface_element_list) )
-     end do
-  else
-     ! from_fields are on the full mesh and we only extract its values
-     ! at the specified surface_elements
-     
-     from_surface_fields=.false.     
-  end if     
-     
-  surface_rhs=0
-  ! Compute contribution of exterior surface integral (boundary condition) to rhs
-  !-----------------------
-  do i=1, size(surface_element_list)
-     f=surface_element_list(i)
-     call transform_facet_to_physical(positions, f, detwei_f)
-     face_mat=-shape_shape(face_shp, face_shp, detwei_f)
-     do k=1, size(from_fields)
-        if (from_surface_fields) then
-           ! we need to use ele_val, where i is the element number in the surface_mesh
-           surface_rhs(:,k,f)=surface_rhs(:,k,f)+ &
-              matmul(face_mat, ele_val(from_fields(k), i))
-        else
-           ! we can simply use face_val with face number f
-           surface_rhs(:,k,f)=surface_rhs(:,k,f)+ &
-              matmul(face_mat, face_val(from_fields(k), f))
-        end if
-     end do
-  end do
-  
-  ! Solution loop
-  !-----------------------
-  if (optimal_ordering) then
-    noit=1
-  else
-    noit=10
-  end if
-  
-  do it=1, noit
-     do i=1, element_count(positions)
-       
-        elm=ordered_elements(i)
-        
-        ! construct diagonal matrix block for this element
-        call transform_to_physical(positions, elm, &
-           shape=ele_shp, dshape=dele_shp, detwei=detwei)
-        ele_mat=shape_vector_dot_dshape(ele_shp, &
-           ele_val_at_quad(vertical_normal,elm), &
-           dele_shp, detwei)
-
-        ! initialise rhs
-        if (present(rhs)) then
-           do k=1, size(to_fields)
-              ele_rhs(:,k)=shape_rhs(ele_shp, detwei*ele_val_at_quad(rhs(k), elm))
-           end do
-        else
-           ele_rhs=0.0
-        end if
-
-        
-        ! then add contribution of surface integrals of incoming
-        ! faces to the rhs and matrix
-        neigh => row_m_ptr(face_normal_gravity, elm)
-        inn => row_val_ptr(face_normal_gravity, elm)
-        faces => ele_faces(positions, elm)
-        do j=1, size(neigh)
-          if (inn(j)<-VERTICAL_INTEGRATION_EPS) then
-            call transform_facet_to_physical(positions, faces(j), &
-                 detwei_f)
-            face_mat=-shape_shape(face_shp, face_shp, detwei_f)*inn(j)
-            
-            face_nds=face_global_nodes(to_fields(1), faces(j))
-            face_lnds => face_local_nodes(to_fields(1)%mesh, faces(j))
-            ele_mat(face_lnds,face_lnds)=ele_mat(face_lnds,face_lnds)+face_mat
-            
-            if (neigh(j)>0) then
-               ! face of element neigh(j), facing elm:
-               f2=ele_face(positions, neigh(j), elm)
-               face_nds2=face_global_nodes(to_fields(1), f2)
-               
-               do k=1, size(to_fields)
-                  ele_rhs(face_lnds,k)=ele_rhs(face_lnds,k)+ &
-                     matmul(face_mat, node_val(to_fields(k), face_nds2))
-               end do
-            else
-               ! note that we've already multiplied with face_mat above, but not with inn(j)
-               ele_rhs(face_lnds,:)=ele_rhs(face_lnds,:)+surface_rhs(:,:,faces(j))*inn(j)
-            end if
-          end if
-        end do
-        
-        call invert(ele_mat)
-        
-        ! compute values for the to_fields:
-        ele_nds => ele_nodes(to_fields(1), elm)
-        do k=1, size(to_fields)
-           call set( to_fields(k), ele_nds, matmul(ele_mat, ele_rhs(:,k)) )
-        end do
-     end do
-  end do
-  
-  call deallocate(face_normal_gravity)
-  
-end subroutine vertical_integration_multiple
-  
-subroutine vertical_extrapolation_module_check_options
-  
-  if (have_option("/geometry/ocean_boundaries")) then
-    if (.not. have_option("/physical_parameters/gravity")) then
-      ewrite(-1,*) "If you select /geometry/ocean_boundaries, you also need to "//&
-        &"set /physical_parameters/gravity"
-      FLExit("Missing gravity!")
-    end if
-  end if
-  
-end subroutine vertical_extrapolation_module_check_options
-
-end module vertical_extrapolation_module

=== modified file 'main/Fluids.F90'
--- main/Fluids.F90	2012-09-03 14:31:15 +0000
+++ main/Fluids.F90	2012-09-04 15:00:37 +0000
@@ -84,7 +84,7 @@
   use discrete_properties_module
   use gls
   use k_epsilon
-  use iceshelf_meltrate_surf_normal
+  use ice_melt_interface
   use halos
   use memory_diagnostics
   use free_surface_module
@@ -406,15 +406,9 @@
 
     call initialise_diagnostics(filename, state)
 
-    ! Initialise ice_meltrate, read constatns, allocate surface, and calculate melt rate
+    ! Initialise melt interface paramerisation
     if (have_option("/ocean_forcing/iceshelf_meltrate/Holland08")) then
-        call melt_surf_init(state(1))
-        call melt_allocate_surface(state(1))
-        call melt_surf_calc(state(1))
-         !BC for ice melt
-          if (have_option('/ocean_forcing/iceshelf_meltrate/Holland08/calculate_boundaries')) then
-            call melt_bc(state(1))
-          endif
+        call melt_interface_initialisation(state(1))
     end if
     
     ! Checkpoint at start
@@ -636,12 +630,13 @@
                 end if
              end if
 
-             ! Calculate the meltrate
-             if(have_option("/ocean_forcing/iceshelf_meltrate/Holland08/") ) then
-                if( (trim(field_name_list(it))=="MeltRate")) then
-                   call melt_surf_calc(state(1))
+            ! Calculate the (ice) interface meltrate using the parameterisation described in Holland et al. 2008 doi:10.1175/2007JCLI1909.1
+            if (have_option("/ocean_forcing/iceshelf_meltrate/Holland08")) then
+                if ((trim(field_name_list(it)) == "MeltRate")) then
+                    call melt_interface_calculate(state(1))
                 endif
-             end if
+            end if
+
 
              call get_option(trim(field_optionpath_list(it))//&
                   '/prognostic/equation[0]/name', &
@@ -713,9 +708,9 @@
              call keps_eddyvisc(state(i))
           end do
           
-          !BC for ice melt
+          ! Boundary conditions for the (ice) melt interface parameterisation
           if (have_option('/ocean_forcing/iceshelf_meltrate/Holland08/calculate_boundaries')) then
-            call melt_bc(state(1))
+            call melt_interface_boundary_condition(state(1))
           endif
           
           if(option_count("/material_phase/scalar_field/prognostic/spatial_discretisation/coupled_cv")>0) then

=== modified file 'main/Makefile.dependencies'
--- main/Makefile.dependencies	2012-03-02 00:25:33 +0000
+++ main/Makefile.dependencies	2012-09-04 15:00:37 +0000
@@ -45,7 +45,7 @@
    ../include/fldebug.mod ../include/foam_flow_module.mod \
    ../include/free_surface_module.mod ../include/global_parameters.mod \
    ../include/gls.mod ../include/goals.mod ../include/halos.mod \
-   ../include/iceshelf_meltrate_surf_normal.mod ../include/implicit_solids.mod \
+   ../include/ice_melt_interface.mod ../include/implicit_solids.mod \
    ../include/k_epsilon.mod ../include/memory_diagnostics.mod \
    ../include/meshdiagnostics.mod ../include/meshmovement.mod \
    ../include/momentum_equation.mod ../include/multimaterial_module.mod \

=== modified file 'manual/embedded_models.tex'
--- manual/embedded_models.tex	2011-11-24 14:46:24 +0000
+++ manual/embedded_models.tex	2012-09-04 15:00:37 +0000
@@ -361,13 +361,6 @@
 therefore requires a single boundary condition: the amount of light incident
 on the water surface. This is a Dirichlet condition on $I$.
 
-As the PAR equation is relatively trivial to solve, the following options
-are recommended:
-\begin{itemize}
-\item Discontinuous discretisation
-\item Solver: gmres iterative method, with no preconditioner.
-\end{itemize}
-
 \subsection{Detritus falling velocity}\label{sec:detritus}
 
 Phytoplankton, zooplankton and nutrient are assumed to be neutrally buoyant

=== modified file 'manual/model_equations.tex'
--- manual/model_equations.tex	2012-08-31 13:12:30 +0000
+++ manual/model_equations.tex	2012-09-04 15:00:37 +0000
@@ -10,7 +10,7 @@
 problem & underlying equations & boundary conditions & other considerations\\
 \hline
 scalar advection & \ref{sec:MP-AdvdifEqn-eqns}, \ref{sec:MP-AdvdifEqn-adv} & \ref{sec:bc_scalar_dirichlet}, \ref{sec:bc_scalar_neumann} & × \\
-scalar advection-diffusion & \ref{sec:MP-AdvdifEqn-eqns}, \ref{sec:MP-AdvdifEqn-diff} & \ref{sec:bc_scalar_dirichlet}, \ref{sec:bc_scalar_neumann}, \ref{sec:bc_scalar_robin} × & ×\\
+scalar advection-diffusion & \ref{sec:MP-AdvdifEqn-eqns}, \ref{sec:MP-AdvdifEqn-diff} & × & ×\\
 momentum equation & \ref{sec:MP-MomEqn} & \ref{sec:bc_vector_dirichlet}, \ref{sec:bc_vector_stress}, \ref{sec:bc_vector_traction}, \ref{sec:FS} & \ref{sec:eqn_extensions} \\
 \hline
 \end{tabular}
@@ -119,21 +119,6 @@
 and substituting in this surface integral to the discretised equation. Note that $q$ is often termed a flux.
 
 
-\subsubsection{Robin condition for a scalar field}\label{sec:bc_scalar_robin}
-\index{boundary conditions!Neumann}
-
-Taking the weak form (applying Green's theorem to the diffusion term) of the advection-diffusion equation \eqref{eq:general_scalar_eqn} leads to a surface integral of the form
-\begin{equation*}
-\int_{\partial\Omega} \phi (\kaptens\nabla c)\cdot\bmn \;d\Gamma,
-\end{equation*}
-where $\phi$ is a test function (see section \ref{chap:numerical_discretisation}).
-The Robin condition is specified by relating the surface diffusive flux $(\kaptens\nabla c)\cdot\bmn$ to an ambient field value $c_{a}$, with an associated surface transfer coefficient $h$. This is given by the relationship
-\begin{equation*}
--(\kaptens\nabla c)\cdot\bmn = h (c - c_{a}),\quad \textrm{on}\quad \partial\Omega,
-\end{equation*}
-where in general $h$ and $c_{a}$ can vary spatially and temporally. This is substituted into the discretised equation to give a surface absorption term given by $hc$ and a surface source given by $-hc_{a}$.
-
-
 \section{Fluid equations}\label{sec:MP-MomEqn}
 
 \index{conservation!equation}
@@ -727,7 +712,7 @@
 
 \subsection{Multi-material simulations}
 \index{multi-material flow}
-The ability to differentiate between regions with distinct material properties is of fundamental importance in the modelling of many physical systems. Two different approaches exist for achieving this: the multi-material approach and the multiphase approach. The multi-material approach is implemented within \fluidity, and the multiphase approach (discussed in the next subsection) is currently under development.
+The ability to differentiate between regions with distinct material properties is of fundamental importance in the modelling of many physical systems. Two different approaches exist for achieving this: the multi-material approach and the multi-phase approach. The multi-material approach is implemented within \fluidity, and the multi-phase approach (discussed in the next subsection) is currently under development.
 In situations where the model can resolve physical mixing of immiscible materials, or where there is no mixing, only one velocity field (and hence one momentum equation) is required to describe the flow of all materials. The \emph{multi-material} approach, considers all materials to be immiscible materials separated by a sharp interface.
 
 In a multi-material approach, the various forms of the conservation equations can be solved for multiple material flows if the point-wise mass density $\rho$ in the equations is defined as the bulk density at each point.  If the flow comprises $n$ materials and the volume fraction of the $i^{th}$ material is denoted $\phi_i$ then the bulk density is given by:
@@ -749,17 +734,14 @@
 \end{equation}
 where the volume fraction field at time zero must be specified.
 
-\subsection{Multiphase simulations}
+\subsection{Multi-phase simulations}
 \label{sec:multiphase_equations}
-\index{multiphase flow}
-Multiphase flows are defined by \cite{prosperettiEtAl2007} to be flows in which two or more phases of matter (solid, liquid, gas, etc) are simultaneously present and are allowed to inter-penetrate. Simple examples include the flow of a fizzy drink which is composed of a liquid and a finite number of gas bubbles, the transportation of solid sediment particles in a river, and the flow of blood cells around the human body.
+\index{multi-phase flow}
+Multi-phase flows are defined by \cite{prosperettiEtAl2007} to be flows in which two or more phases of matter (solid, liquid, gas, etc) are simultaneously present and are allowed to inter-penetrate. Simple examples include the flow of a fizzy drink which is composed of a liquid and a finite number of gas bubbles, the transportation of solid sediment particles in a river, and the flow of blood cells around the human body.
 
 Further to the above definition, each phase is classed as either \textit{continuous} or \textit{dispersed}, where a continuous phase is a connected liquid or gas substance in which dispersed phases (comprising a finite number of solid particles, liquid droplets and/or gas bubbles) may be immersed \citep{croweEtAl1998}.
 
-To enable the mixing and inter-penetration of phases, a separate velocity field (and hence a separate momentum equation) is assigned to each one and solved for. Extra terms are then included to account for inter-phase interactions. Furthermore, the model currently assumes no mass transfer between phases, incompressible flow, and a common pressure field $p$ so that only one continuity equation is used. The form of the governing equations depends on whether we are dealing with incompressible or compressible flow.
-
-\subsubsection{Incompressible flow}
-For an incompressible multiphase flow, the continuity equation and momentum equation for phase $i$ (based on the derivation in \cite{ishii1975}, written in non-conservative form) are:
+To enable the mixing and inter-penetration of phases, a separate velocity field (and hence a separate momentum equation) is assigned to each one and solved for. Extra terms are then included to account for inter-phase interactions. Furthermore, the model currently assumes no mass transfer between phases, incompressible flow, and a common pressure field $p$ so that only one continuity equation is used. Thus, the continuity equation and momentum equation for phase $i$ (based on the derivation in \cite{ishii1975}, written in non-conservative form) are:
 \begin{equation}
 \sum_{i=1}^N\nabla\cdot\left(\alpha_i\mathbf{u}_i\right) = 0,
 \end{equation}
@@ -772,19 +754,6 @@
 \end{equation}
 where $\mathbf{u}_i$, $\rho_i$, $\tautens_i$ and $\alpha_i$ are the velocity, density, viscous stress tensor and volume fraction of phase $i$ respectively, and $\mathbf{F}_i$ represents the forces imposed on phase $i$ by the other $N-1$ phases. Details of momentum transfer terms are given in Chapter \ref{chap:configuration}.
 
-\subsubsection{Compressible flow}
-Fluidity currently only supports the case where the continuous phase is compressible. All other phases (the dispersed phases) are incompressible. The momentum equation is the same as the one given above, but the continuity equation becomes:
-\begin{equation}
-\alpha_1\frac{\partial\rho_1}{\partial t} + \nabla\cdot\left(\alpha_1\rho_1\mathbf{u}_1\right) + \rho_1\left(\sum_{i=2}^N\nabla\cdot\left(\alpha_i\mathbf{u}_i\right)\right) = 0,
-\end{equation}
-where the compressible phase has an index of $i=1$ here.
-
-An internal energy equation may also need to be solved depending on the equation of state used for the compressible phase. This is given by:
-\begin{equation}
-\alpha_i\rho_i\frac{\partial e_i}{\partial t} + \alpha_i\rho_i\mathbf{u}_i\cdot\nabla e_i = -p\nabla\cdot\left(\alpha_i\mathbf{u}_i\right) + Q_i
-\end{equation}
-where $e_i$ is the internal energy and $Q_i$ is an inter-phase energy transfer term described in Chapter \ref{chap:configuration}.
-
 \subsection{Porous Media Darcy Flow}
 \label{sec:porous_media_darcy_flow_equations}
 \index{porous media darcy flow}

=== modified file 'manual/numerical_discretisation.tex'
--- manual/numerical_discretisation.tex	2012-07-26 15:38:42 +0000
+++ manual/numerical_discretisation.tex	2012-09-04 15:00:37 +0000
@@ -9,12 +9,11 @@
 $\Omega$ is defined as $\dOmega$
 %(NB. $\Gamma$ is another common symbol for a boundary) 
 and can be split into different sections, e.g.
-$\partial\Omega = \dOmega^{N}\cup \dOmega^{R}\cup \dOmega^{D}$
+$\partial\Omega = \dOmega^{N}\cup \dOmega^{D}$
 where $\dOmega^{N}$ is that part of the boundary over
-which Neumann conditions are applied, $\dOmega^{R}$ is the part of the 
-boundary over which Robin conditions are applied and $\dOmega^{D}$ is that part of the boundary over
+which Neumann conditions are applied and $\dOmega^{D}$ is that part of the boundary over
 which Dirichlet conditions are applied. For clarity a subscript showing the field in question will
-be given on the $N$, $R$ and $D$ characters.
+be given on the $N$ and $D$ characters.
 
 The unit vector \vec{n} is always assumed to be the
 outward facing normal vector to the domain.
@@ -31,11 +30,10 @@
 \end{equation}
 
 We now define a partition of the boundary such that
-$\partial\Omega = \dOmega^{N}\cup \dOmega^{R}\cup \dOmega^{D}$, and impose boundary conditions on $c$:
+$\partial\Omega = \dOmega^{N}\cup \dOmega^{D}$, and impose boundary conditions on $c$:
 \begin{gather}
   \label{eq:cboundary}
   \vec{n}\cdot\tensor{\kappa}\cdot\grad c=g_{N_c} \quad\textrm{on}\quad \dOmega^{N_c},\\
-  -\vec{n}\cdot\tensor{\kappa}\cdot\grad c=h(c-c_{a})\quad\textrm{on}\quad \dOmega^{R_c},\\
   c=g_{D_c} \quad\textrm{on}\quad \dOmega^{D_c}.
 \end{gather}
 
@@ -371,7 +369,7 @@
 
 The SUPG implementation in \fluidity\ does not modify the test function derivatives
 or the face test functions. Hence the SUPG implementation is only consistent
-for degree one elements with no non-zero Neumann, Robin or weak Dirichlet boundary
+for degree one elements with no non-zero Neumann or weak Dirichlet boundary
 conditions.
 
 SUPG is considered ready for production use for scalar advection-diffusion
@@ -435,25 +433,6 @@
 boundary term \eqref{eq:missing_boundary_term} with an integral of $\phi~g_N$
 over the boundary.
 
-The Robin boundary condition
-\begin{equation*}
-  -\vec{n}\cdot\tensor{\kappa}\cdot\grad c=h(c-c_{a}),
-\end{equation*}
-where $h$ and $c_{a}$ can be any prescribed function on the boundary $\partial\Omega$, is also imposed 
-weakly. As for the Neumann boundary condition, weighting the Robin boundary condition with a test function 
-and integrating over the relevant domain boundary gives the weak form:
-\begin{equation}\label{eq:weak_robin}
-  - \int_{\partial\Omega} \phi~\vec{n}\cdot\tensor{\kappa}\cdot\grad c=
-    \int_{\partial\Omega} \phi~h(c-c_{a}), \quad \text{for all }\phi.
-\end{equation}
-Thus \eqref{eq:weak_robin} can be used to replace the missing
-boundary term \eqref{eq:missing_boundary_term} with an integral of $\phi~h(c-c_{a})$
-over the relevant part of the boundary. The two terms within the integral can be treated slightly differently. 
-The term $\phi~hc_{a}$ is always included in the right hand side of the linear system. 
-The term $\phi~hc$ is effectively a surface absorption term where any implicit contributions 
-are included in the left had side matrix, while any explicit contributions are included into the 
-right hand side of the linear system.
-
 \index{boundary conditions!Dirichlet!weakly imposed}
 In a similar way, a weakly imposed Dirichlet boundary condition can be related to an
 integration by parts of the advection term. Let us consider a pure advection problem

=== modified file 'manual/parameterisations.tex'
--- manual/parameterisations.tex	2012-08-30 14:23:11 +0000
+++ manual/parameterisations.tex	2012-09-04 15:00:37 +0000
@@ -83,18 +83,6 @@
 \end{equation}
 where $c_\mu^0$ is a model constant used to make $\Psi$ identifiable with any of the transitional
 models, e.g. $kl$, $\epsilon$, and $\omega$, by adopting the values shown in Table \ref{tab:glsparams} \citep{umlauf2003}.
-
-There is also the option to add an extra term to account for various
-oceanic parameters, such an internal waves breaking. This is based on the NEMO
-ocean model and takes a user-defined percentage of the surface $k$ and adds it down-depth
-using an exponential profile:
-\begin{equation}
-    k_z = {k_z}_o + p*k_{\mathrm{sur}} * \exp{\left( -z / l_k \right)}
-\end{equation}
-where $k_z$ is the new turbulent kinetic energy value at depth, $z$, ${k_z}_o$ is the original TKE, $k_{\mathrm{sur}}$
-is the surface TKE, $p$ is a constant, and $l_k$ is a lengthscale. The options for this can be found in
-\option{\ldots/subgridscale\_parameterisations/gls/ocean\_parameterisation}.
-
 The second equation is:
 \begin{equation}
 \frac{\partial \Psi}{\partial t} + \bmu_i\frac{\partial \Psi}{\partial x_i} =

=== modified file 'manual/visualisation_and_diagnostics.tex'
--- manual/visualisation_and_diagnostics.tex	2012-08-09 22:08:12 +0000
+++ manual/visualisation_and_diagnostics.tex	2012-09-04 15:00:37 +0000
@@ -66,7 +66,6 @@
 Limitations: Requires the diagnostic \option{scalar\_field::IsopycnalCoordinate} and (therefore) is not parallelised.
 \item[BulkMaterialPressure:]Calculates the bulk material pressure based on the MaterialDensity and MaterialVolumeFraction (and MaterialInternalEnergy if appropriate) for the equation of state of all materials.
 \item[CFLNumber:]CFLNumber as defined on the co-ordinate mesh. It is calculated as $\triangle t \vec{u} \mat{J}^{-1}$ where $\triangle t$ is the timestep, $\vec{u}$ the velocity and $\mat{J}$ the Jacobian. 
-\item[CompressibleContinuityResidual:]Computes the residual of the compressible multiphase continuity equation. Used only in compressible multiphase flow simulations. 
 \item[ControlVolumeCFLNumber:]CFL Number as defined on a control volume mesh. It is calculated as $\triangle t \frac{1}{V} \int _{cv_{face}} \vec{u}$ where $\triangle t$ is the timestep and $\vec{u}$ the velocity. The integral is taken over the faces of the control volume and $V$ is the volume of the control volume.
 \item[ControlVolumeDivergence:]Divergence of the velocity field where the divergence operator is defined using the control volume $\mathrm{C}^\mathrm{T}$ matrix. This assumes that the test space is discontinuous control volumes.
 \item[CVMaterialDensityCFLNumber:]Courant Number as defined on a control volume mesh and incorporating the MaterialDensity. Requires a MaterialDensity field!

=== modified file 'parameterisation/Gls_vertical_turbulence_model.F90'
--- parameterisation/Gls_vertical_turbulence_model.F90	2012-06-08 20:26:43 +0000
+++ parameterisation/Gls_vertical_turbulence_model.F90	2012-09-04 15:00:37 +0000
@@ -42,7 +42,6 @@
   use Coordinates
   use FLDebug
   use fefields
-  use vertical_extrapolation_module
 
   implicit none
 
@@ -355,13 +354,10 @@
     type(tensor_field), pointer      :: tke_diff, background_diff
     type(scalar_field), pointer      :: tke
     real                             :: prod, buoyan, diss
-    integer                          :: i, stat, ele
+    integer                          :: i, stat
     character(len=FIELD_NAME_LEN)    :: bc_type
     type(scalar_field), pointer      :: scalar_surface
     type(vector_field), pointer      :: positions
-    type(scalar_field), pointer      :: lumped_mass
-    type(scalar_field)               :: inverse_lumped_mass
-
 
     ! Temporary tensor to hold  rotated values if on the sphere (note: must be a 3x3 mat)
     real, dimension(3,3) :: K_M_sphere_node
@@ -371,10 +367,6 @@
     ! Get N^2 and M^2 -> NN2 and MM2
     call gls_buoyancy(state)
 
-    do i=1,NNodes_sur
-        call set(NN2,top_surface_nodes(i),0.0)
-    end do
-
     ! calculate stability function
     call gls_stability_function(state)
 
@@ -384,35 +376,26 @@
     tke_diff => extract_tensor_field(state, "GLSTurbulentKineticEnergyDiffusivity")
     positions => extract_vector_field(state, "Coordinate")
 
-    ! Create a local_tke in which we can mess about with the surface values
-    ! for creating some of the diagnostic values later
+    ! swap state tke for our local one - the state tke has had the upper and
+    ! lower boundaries ammended with the Dirichlet condition as in Warner et al
+    ! 2005. The local TKE is the non-amended version. This ensures a
+    ! non-ill-posed problem.
     call set(tke,local_tke)
-    ! Assembly loop
-    call zero(P)
-    call zero(B)
-    call allocate(inverse_lumped_mass, P%mesh, "InverseLumpedMass")
-    lumped_mass => get_lumped_mass(state, P%mesh)
-    call invert(lumped_mass, inverse_lumped_mass)
-    ! create production terms, P, B
-    do ele=1, ele_count(P)
-        call assemble_tke_prodcution_terms(ele, P, B, mesh_dim(P))
-    end do
-    call scale(P,inverse_lumped_mass)
-    call scale(B,inverse_lumped_mass)
-    call deallocate(inverse_lumped_mass)
 
-    call zero(source)
-    call zero(absorption)
-    do ele = 1, ele_count(tke)
-        call assemble_kk_src_abs(ele,tke, mesh_dim(tke))
-    end do
-    call allocate(inverse_lumped_mass, tke%mesh, "InverseLumpedMass")
-    lumped_mass => get_lumped_mass(state, tke%mesh)
-    call invert(lumped_mass, inverse_lumped_mass)
-    ! source and absorption terms are set, apart from the / by lumped mass
-    call scale(source,inverse_lumped_mass)
-    call scale(absorption,inverse_lumped_mass)
-    call deallocate(inverse_lumped_mass)
+    do i=1,nNodes
+        prod = node_val(K_M,i)*node_val(MM2,i) 
+        buoyan = -1.*node_val(K_H,i)*node_val(NN2,i)
+        diss = node_val(eps,i)
+        if (prod+buoyan.gt.0) then
+            call set(source,i, prod + buoyan)
+            call set(absorption,i, diss / node_val(tke,i))
+        else
+            call set(source,i, prod)
+            call set(absorption,i, (diss-buoyan)/node_val(tke,i))
+        end if
+        call set(P, i, prod)
+        call set(B, i, buoyan)
+    enddo
 
     ! set diffusivity for tke
     call zero(tke_diff)
@@ -457,66 +440,6 @@
        call set(scalarField,absorption)  
     end if
 
-    contains
-
-    subroutine assemble_tke_prodcution_terms(ele, P, B, dim)
-
-        integer, intent(in)               :: ele, dim
-        type(scalar_field), intent(inout) :: P, B
-
-        real, dimension(ele_loc(P,ele),ele_ngi(P,ele),dim)  :: dshape_P
-        real, dimension(ele_ngi(P,ele))                     :: detwei
-        real, dimension(ele_loc(P,ele))                     :: rhs_addto_vel, rhs_addto_buoy
-        type(element_type), pointer                         :: shape_p
-        integer, pointer, dimension(:)                      :: nodes_p
-    
-        nodes_p => ele_nodes(p, ele)
-        shape_p => ele_shape(p, ele)
-        call transform_to_physical( positions, ele, shape_p, dshape=dshape_p, detwei=detwei )
-
-        ! Shear production term:
-        rhs_addto_vel = shape_rhs(shape_p, detwei*ele_val_at_quad(K_M,ele)*ele_val_at_quad(MM2,ele))
-        ! Buoyancy production term:
-        rhs_addto_buoy = shape_rhs(shape_p, -detwei*ele_val_at_quad(K_H,ele)*ele_val_at_quad(NN2,ele))
- 
-        call addto(P, nodes_p, rhs_addto_vel)
-        call addto(B, nodes_p, rhs_addto_buoy)
-    
-    end subroutine assemble_tke_prodcution_terms
-
-    subroutine assemble_kk_src_abs(ele, kk, dim)
-
-        integer, intent(in)            :: ele, dim
-        type(scalar_field), intent(in) :: kk
-
-        real, dimension(ele_loc(kk,ele),ele_ngi(kk,ele),dim) :: dshape_kk
-        real, dimension(ele_ngi(kk,ele))                     :: detwei
-        real, dimension(ele_loc(kk,ele))                     :: rhs_addto_disip, rhs_addto_src
-        type(element_type), pointer                          :: shape_kk
-        integer, pointer, dimension(:)                       :: nodes_kk
-    
-        nodes_kk => ele_nodes(kk, ele)
-        shape_kk => ele_shape(kk, ele)
-        call transform_to_physical( positions, ele, shape_kk, dshape=dshape_kk, detwei=detwei )
-
-        ! ROMS and GOTM hide the absorption term in the source if the
-        ! total is > 0. Kinda hard to do that over an element (*what* should
-        ! be > 0?). So we don't pull this trick. Doesn't seem to make a
-        ! difference to anything.
-        rhs_addto_src = shape_rhs(shape_kk, detwei * (&
-                                 (ele_val_at_quad(P,ele))) &
-                                 )
-        rhs_addto_disip = shape_rhs(shape_kk, detwei * ( &
-                                 (ele_val_at_quad(eps,ele) - &
-                                 ele_val_at_quad(B,ele)) / &
-                                 ele_val_at_quad(tke,ele)) &
-                                 )
-           
-        call addto(source, nodes_kk, rhs_addto_src)
-        call addto(absorption, nodes_kk, rhs_addto_disip)
-
-
-    end subroutine assemble_kk_src_abs
 
 end subroutine gls_tke
 
@@ -532,19 +455,9 @@
     type(tensor_field), pointer      :: psi_diff, background_diff
     real                             :: prod, buoyan,diss,PsiOverTke
     character(len=FIELD_NAME_LEN)    :: bc_type
-    integer                          :: i, stat, ele
+    integer                          :: i, stat
     type(scalar_field), pointer      :: scalar_surface
     type(vector_field), pointer      :: positions
-    type(scalar_field), pointer      :: lumped_mass
-    type(scalar_field)               :: inverse_lumped_mass, vel_prod, buoy_prod
-    ! variables for the ocean parameterisation
-    type(csr_matrix)                   :: face_normal_gravity
-    integer, dimension(:), allocatable :: ordered_elements
-    logical, dimension(:), allocatable :: node_list
-    logical                            :: got_surface
-    real                               :: lengthscale, percentage, tke_surface
-    type(scalar_field), pointer        :: distanceToTop, distanceToBottom
-    type(vector_field), pointer        :: vertical_normal
 
     ! Temporary tensor to hold  rotated values (note: must be a 3x3 mat)
     real, dimension(3,3)             :: psi_sphere_node
@@ -558,16 +471,26 @@
     psi => extract_scalar_field(state, "GLSGenericSecondQuantity")
     psi_diff => extract_tensor_field(state, "GLSGenericSecondQuantityDiffusivity")
     positions => extract_vector_field(state, "Coordinate")
-    vertical_normal => extract_vector_field(state, "GravityDirection")
-    call allocate(vel_prod, psi%mesh, "_vel_prod_psi")
-    call allocate(buoy_prod, psi%mesh, "_buoy_prod_psi")
+
     ewrite(2,*) "In gls_psi: setting up"
 
     ! store the tke in an internal field and then
     ! add the dirichlet conditions to the upper and lower surfaces. This
     ! helps stabilise the diffusivity (e.g. rapid heating cooling of the surface
-    ! can destabilise the run)
+    ! can destabilise the run) but also is then done for output so we match
+    ! other models which also play this trick, c.f. Warner et al 2005.
     call set(local_tke,tke)
+    ! Call the bc code, but specify we want dirichlet
+    call gls_tke_bc(state, 'dirichlet')
+    ! copy the values onto the mesh using the global node id
+    do i=1,NNodes_sur
+        call set(tke,top_surface_nodes(i),node_val(top_surface_values,i))
+        call set(tke_old,top_surface_nodes(i),node_val(top_surface_values,i))   
+    end do   
+    do i=1,NNodes_bot
+        call set(tke,bottom_surface_nodes(i),node_val(bottom_surface_values,i))
+        call set(tke_old,bottom_surface_nodes(i),node_val(bottom_surface_values,i))
+    end do   
 
     ! clip at k_min
     do i=1,nNodes
@@ -576,76 +499,41 @@
         call set(local_tke,i, max(node_val(local_tke,i),k_min))
     end do
 
-    ! This is the extra term meant to add in internal wave breaking and the
-    ! like. Based on a similar term in NEMO
-    if (have_option("/material_phase[0]/subgridscale_parameterisations/GLS/ocean_parameterisation")) then
-
-        distanceToTop => extract_scalar_field(state, "DistanceToTop")
-        distanceToBottom => extract_scalar_field(state, "DistanceToBottom")
-
-        call get_option("/material_phase[0]/subgridscale_parameterisations/GLS/ocean_parameterisation/lengthscale",lengthscale)
-        call get_option("/material_phase[0]/subgridscale_parameterisations/GLS/ocean_parameterisation/percentage",percentage)
-        ewrite(2,*) "Computing extra ocean parameterisation"
-        allocate(node_list(NNodes))
-        node_list = .false.
-        ! create gravity face normal
-        call compute_face_normal_gravity(face_normal_gravity, positions, vertical_normal)
-        allocate(ordered_elements(size(face_normal_gravity,1)))
-
-        ! create an element ordering from the mesh that moves vertically downwards
-        call vertical_element_ordering(ordered_elements, face_normal_gravity)
-        ! I assume the above fails gracefully if the mesh isn't suitable, hence
-        ! no options checks are carried out
-
-        got_surface = .false.
-        do i=1, size(ordered_elements)
-            if (.not. got_surface) then
-                ! First time around grab, the surface TKE
-                ! for this column
-                tke_surface = maxval(ele_val(tke,i))
-                got_surface = .true.
-            end if
-            if (got_surface) then
-                call ocean_tke(i,tke,distanceToTop,lengthscale,percentage,tke_surface, node_list)
-            end if
-            if (minval(ele_val(distanceToBottom,i)) < 1e-6) then
-                got_surface = .false.
-            end if
-        end do
-
-        deallocate(ordered_elements)
-        call deallocate(face_normal_gravity)
-        deallocate(node_list)
-
-    end if
-
-    call allocate(inverse_lumped_mass, psi%mesh, "InverseLumpedMass")
-    lumped_mass => get_lumped_mass(state, psi%mesh)
-    call invert(lumped_mass, inverse_lumped_mass)
-    ! Set Psi from previous timestep
+    ! re-construct psi at "old" timestep
     do i=1,nNodes
         call set(psi,i, cm0**gls_p * node_val(tke_old,i)**gls_m * node_val(ll,i)**gls_n)
         call set(psi,i,max(node_val(psi,i),psi_min))
     end do
 
     ewrite(2,*) "In gls_psi: computing RHS"
-    call zero(vel_prod)
-    call zero(buoy_prod)
-    call zero(source)
-    call zero(absorption)    
-    do ele = 1, ele_count(psi)
-        call assemble_production_terms_psi(ele, vel_prod, buoy_prod, psi, mesh_dim(psi))
-    end do
-    call scale(vel_prod,inverse_lumped_mass)
-    call scale(buoy_prod,inverse_lumped_mass)
-
-    do ele = 1, ele_count(psi)
-        call assemble_psi_src_abs(ele, psi, tke_old, mesh_dim(psi))
-    end do
-    call scale(source,inverse_lumped_mass)
-    call scale(absorption,inverse_lumped_mass)
-    call deallocate(inverse_lumped_mass)
-    
+    ! compute RHS
+    do i=1,nNodes
+
+        ! compute production terms in psi-equation
+        if (node_val(B,i).gt.0) then ! note that we have already set B when setting up the RHS for TKE
+            cPsi3=cPsi3_plus  ! unstable strat
+        else
+            cPsi3=cPsi3_minus ! stable strat
+        end if
+
+        ! compute production terms in psi-equation
+        PsiOverTke  = node_val(psi,i)/node_val(tke_old,i)
+        prod        = cPsi1*PsiOverTke*node_val(P,i)
+        buoyan      = cPsi3*PsiOverTke*node_val(B,i)
+        diss        = cPsi2*PsiOverTke*node_val(eps,i)*node_val(Fwall,i)
+        if (prod+buoyan.gt.0) then
+            src = prod+buoyan
+            absn = diss/node_val(psi,i)
+            call set(source,i, src)
+            call set(absorption,i, absn)
+        else
+            src = prod
+            absn = (diss-buoyan)/node_val(psi,i)
+            call set(source,i, src)
+            call set(absorption,i, absn)
+        end if
+    end do
+
     ewrite(2,*) "In gls_psi: setting diffusivity"
     ! Set diffusivity for Psi
     call zero(psi_diff)
@@ -689,150 +577,6 @@
        call set(scalarField,absorption)  
     end if
 
-    call deallocate(vel_prod)
-    call deallocate(buoy_prod)
-
-
-    contains 
-    
-    subroutine ocean_tke(ele, tke, distanceToTop, lengthscale, percentage, tke_surface, nodes_done)
-        type(scalar_field),pointer, intent(in)  :: distanceToTop
-        type(scalar_field),pointer, intent(out) :: tke
-        real, intent(in)                        :: lengthscale, percentage, tke_surface
-        integer, intent(in)                     :: ele
-        logical, dimension(:), intent(inout)    :: nodes_done
-
-        integer, dimension(:), pointer  :: element_nodes
-        integer                         :: i, node
-        real                            :: current_TKE, depth
-
-        element_nodes => ele_nodes(tke, ele)
-
-
-        ! smooth out TKE according to length scale
-        do i = 1, size(element_nodes)
-            node = element_nodes(i)
-            depth = node_val(distanceToTop,node)
-            current_TKE = node_val(tke,node)
-            if (nodes_done(node)) then
-                cycle
-            end if
-            current_TKE = current_TKE + &
-                           & percentage*TKE_surface * EXP( -depth / lengthscale )
-            call set(tke,node,current_TKE)  
-
-            nodes_done(node) = .true.
-        end do
-
-    end subroutine ocean_tke
-
-    subroutine reconstruct_psi(ele, psi, dim)
-
-        integer, intent(in)               :: ele, dim
-        type(scalar_field), intent(inout) :: psi
-
-        real, dimension(ele_loc(psi,ele),ele_ngi(psi,ele),dim) :: dshape_psi
-        real, dimension(ele_ngi(psi,ele))                      :: detwei
-        real, dimension(ele_loc(psi,ele))                      :: rhs_addto
-        type(element_type), pointer                            :: shape_psi
-        integer, pointer, dimension(:)                         :: nodes_psi
-    
-        nodes_psi => ele_nodes(psi, ele)
-        shape_psi => ele_shape(psi, ele)
-        call transform_to_physical( positions, ele, shape_psi, dshape=dshape_psi, detwei=detwei )
-
-        rhs_addto = shape_rhs(shape_psi, detwei* &
-                                         (cm0**gls_p) * &
-                                         ele_val_at_quad(tke_old,ele)**gls_m * &
-                                         ele_val_at_quad(ll,ele)**gls_n)
-
-        call addto(psi, nodes_psi, rhs_addto)
-
-    end subroutine reconstruct_psi
-
-    subroutine assemble_production_terms_psi(ele, vel_prod, buoy_prod, psi, dim)
-
-
-        integer, intent(in)               :: ele, dim
-        type(scalar_field), intent(inout) :: psi, vel_prod, buoy_prod
-
-        real, dimension(ele_loc(psi,ele),ele_ngi(psi,ele),dim) :: dshape_psi
-        real, dimension(ele_ngi(psi,ele))                      :: detwei
-        real, dimension(ele_loc(psi,ele))                      :: rhs_addto_vel, rhs_addto_buoy
-        real, dimension(ele_ngi(psi,ele))                      :: cPsi3
-        type(element_type), pointer                            :: shape_psi
-        integer, pointer, dimension(:)                         :: nodes_psi
-        
-        nodes_psi => ele_nodes(psi, ele)
-        shape_psi => ele_shape(psi, ele)
-        call transform_to_physical( positions, ele, shape_psi, dshape=dshape_psi, detwei=detwei )
-
-        ! Buoyancy production term:
-        ! First we need to work out if cPsi3 is for stable or unstable
-        ! stratification
-        where(ele_val_at_quad(B,ele) .gt. 0.0)
-            cPsi3 = cPsi3_plus  ! unstable strat
-        elsewhere
-            cPsi3 = cPsi3_minus ! stable strat
-        end where
-        rhs_addto_buoy = shape_rhs(shape_psi, detwei*(cPsi3*ele_val_at_quad(B,ele)*&
-                                   (ele_val_at_quad(psi, ele)/ele_val_at_quad(local_tke,ele))))
-
-        ! shear production term:
-        rhs_addto_vel = shape_rhs(shape_psi, detwei*(cPsi1*ele_val_at_quad(P,ele)*&
-                                    (ele_val_at_quad(psi, ele)/ele_val_at_quad(local_tke,ele))))
- 
-        call addto(vel_prod, nodes_psi, rhs_addto_vel)
-        call addto(buoy_prod, nodes_psi, rhs_addto_buoy)
-
-
-    end subroutine assemble_production_terms_psi
-
-    subroutine assemble_psi_src_abs(ele, psi, tke, dim)
-
-        integer, intent(in)            :: ele, dim
-        type(scalar_field), intent(in) :: psi, tke
-
-        real, dimension(ele_loc(psi,ele),ele_ngi(psi,ele),dim) :: dshape_psi
-        real, dimension(ele_ngi(psi,ele))                      :: detwei
-        real, dimension(ele_loc(psi,ele))                      :: rhs_addto_src, rhs_addto_disip
-        type(element_type), pointer                            :: shape_psi
-        integer, pointer, dimension(:)                         :: nodes_psi
-    
-        nodes_psi => ele_nodes(psi, ele)
-        shape_psi => ele_shape(psi, ele)
-        call transform_to_physical( positions, ele, shape_psi, dshape=dshape_psi, detwei=detwei )
-
-        where (ele_val_at_quad(vel_prod,ele) + ele_val_at_quad(buoy_prod,ele) .gt. 0)
-            rhs_addto_src = shape_rhs(shape_psi, detwei* ( &
-                                (ele_val_at_quad(vel_prod,ele)) + &
-                                (ele_val_at_quad(buoy_prod,ele)) &
-                                 ) & !detwei
-                                 ) ! shape_rhs
-            rhs_addto_disip = shape_rhs(shape_psi, detwei* (&
-                                  (cPsi2*ele_val_at_quad(eps,ele) * &
-                                  (ele_val_at_quad(Fwall,ele)*ele_val_at_quad(psi, ele)/ele_val_at_quad(tke,ele))) / &
-                                  ele_val_at_quad(psi,ele) &
-                                  ) & ! detwei
-                                  ) !shape_rhs
-        elsewhere
-            rhs_addto_src = shape_rhs(shape_psi, detwei * ( &
-                                (ele_val_at_quad(vel_prod,ele)) &
-                                 ) & !detwei
-                                 ) !shape_rhs
-            rhs_addto_disip = shape_rhs(shape_psi, detwei * (&
-                                  ((cPsi2*ele_val_at_quad(eps,ele) * &
-                                  (ele_val_at_quad(Fwall,ele)*ele_val_at_quad(psi, ele)/ele_val_at_quad(tke,ele))) - &! disipation term
-                                  (ele_val_at_quad(buoy_prod,ele)))/ & ! buoyancy term
-                                  ele_val_at_quad(psi,ele)&
-                                  ) & !detwei 
-                                  ) !shape_rhs
-        end where
-        call addto(source, nodes_psi, rhs_addto_src)
-        call addto(absorption, nodes_psi, rhs_addto_disip)
-
-    end subroutine assemble_psi_src_abs
-
 end subroutine gls_psi
 
 
@@ -872,15 +616,27 @@
     velocity => extract_vector_field(state, "Velocity")
 
     call allocate(tke, tke_state%mesh, name="MyLocalTKE")
-    !if (gls_n > 0) then
+    if (gls_n > 0) then
         ! set the TKE to use below to the unaltered TKE
-        ! with no changes to the upper/lower surfaces
+        ! with no chnages tothe upper/lower surfaces
         ! Applies to k-kl model only
-    !    call set (tke, local_tke)
-    !else
+        call set (tke, local_tke)
+    else
         ! Use the altered TKE to get the surface diffusivity correct
         call set (tke, tke_state)
-    !end if
+    end if
+
+    ! call the bc code, but specify we want dirichlet
+    if (gls_n < 0) then
+        call gls_psi_bc(state, 'dirichlet')
+        ! copy the values onto the mesh using the global node id
+        do i=1,NNodes_sur
+            call set(psi,top_surface_nodes(i),node_val(top_surface_values,i))
+        end do
+        do i=1,NNodes_bot
+            call set(psi,bottom_surface_nodes(i),node_val(bottom_surface_values,i))
+        end do
+    end if
 
     exp1 = 3.0 + gls_p/gls_n
     exp2 = 1.5 + gls_m/gls_n
@@ -913,12 +669,12 @@
 
         ! compute dissipative scale
         call set(ll,i,cde*sqrt(tke_cur**3.)/node_val(eps,i))
-        !if (gls_n > 0) then
-        !    if (node_val(NN2,i) > 0) then
-        !        limit = sqrt(0.56 * tke_cur / node_val(NN2,i))
-        !        call set(ll,i,min(limit,node_val(ll,i)))
-        !    end if
-        !end if
+        if (gls_n > 0) then
+            if (node_val(NN2,i) > 0) then
+                limit = sqrt(0.56 * tke_cur / node_val(NN2,i))
+                call set(ll,i,min(limit,node_val(ll,i)))
+            end if
+        end if
 
     end do
 
@@ -1684,24 +1440,20 @@
     case("neumann")
         do i=1,NNodes_sur
             ! GOTM Boundary
-            value = -(gls_n*(cm0**(gls_p+1.))*(kappa**(gls_n+1.)))/sigma_psi     &
-                     *node_val(top_surface_tke_values,i)**(gls_m+0.5)*(z0s(i))**gls_n  
-            ! Warner 2005 - left here for posterity and debugging
-            !value = -gls_n*(cm0**(gls_p))*(node_val(top_surface_tke_values,i)**gls_m)* &
-            !         (kappa**gls_n)*(z0s(i)**(gls_n-1))*((node_val(top_surface_km_values,i)/sigma_psi))
+            !value = -(gls_n*(cm0**(gls_p+1.))*(kappa**(gls_n+1.)))/sigma_psi     &
+            !         *node_val(top_surface_values,i)**(gls_m+0.5)*(z0s(i))**gls_n  
+            ! Warner 2005
+            value = -gls_n*(cm0**(gls_p))*(node_val(top_surface_tke_values,i)**gls_m)* &
+                     (kappa**gls_n)*(z0s(i)**(gls_n-1))*((node_val(top_surface_km_values,i)/sigma_psi))
             call set(top_surface_values,i,value)
         end do
         do i=1,NNodes_bot
-            if (u_taub_squared(i) < 1e-16) then
-                value = 0.0
-            else
-                ! GOTM Boundary
-                value = - gls_n*cm0**(gls_p+1.)*(kappa**(gls_n+1.)/sigma_psi)      &
-                           *node_val(bottom_surface_tke_values,i)**(gls_m+0.5)*(z0b(i))**gls_n
-                ! Warner 2005 - as above
-                !value = gls_n*cm0**(gls_p)*node_val(bottom_surface_tke_values,i)**(gls_m)* &
-                !         kappa**gls_n*(z0b(i)**(gls_n-1))*(node_val(bottom_surface_km_values,i)/sigma_psi)
-            end if
+            ! GOTM Boundary
+            !value = - gls_n*cm0**(gls_p+1.)*(kappa**(gls_n+1.)/sigma_psi)      &
+            !           *node_val(bottom_surface_tke_values,i)**(gls_m+0.5)*(z0b(i))**gls_n
+            ! Warner 2005
+            value = gls_n*cm0**(gls_p)*node_val(bottom_surface_tke_values,i)**(gls_m)* &
+                     kappa**gls_n*(z0b(i)**(gls_n-1))*(node_val(bottom_surface_km_values,i)/sigma_psi)
             call set(bottom_surface_values,i,value)
         end do
     case("dirichlet")
@@ -1741,7 +1493,7 @@
     integer                              :: nobcs
     integer                              :: i,ii, MaxIter
     real                                 :: rr
-    real                                 :: charnock_val=18500.
+    real                                 :: charnock_val=1400.
     character(len=OPTION_PATH_LEN)       :: bctype
     type(vector_field), pointer          :: wind_surface_field, positions, velocity
     type(scalar_field), pointer          :: tke
@@ -1835,15 +1587,12 @@
 
         do i=1,NNodes_bot
             temp_vector_3D = node_val(bottom_velocity,i)
-            u_taub = sqrt(temp_vector_3D(1)**2+temp_vector_3D(2)**2+temp_vector_3D(3)**2)
-            if (u_taub <= 1e-12) then
-                z0b(i) = z0s_min
-            else
+            u_taub = max(1e-12,sqrt(temp_vector_3D(1)**2+temp_vector_3D(2)**2+temp_vector_3D(3)**2))
 
 
             !  iterate bottom roughness length MaxIter times
             do ii=1,MaxIter
-                    z0b(i)=(1e-7/max(1e-6,u_taub)+0.03*0.1)
+                z0b(i)=1e-7/max(1e-6,u_taub)+0.03*0.1
 
                 ! compute the factor r
                 rr=kappa/log(z0b(i))
@@ -1852,7 +1601,6 @@
                 u_taub = rr*sqrt(temp_vector_3D(1)**2+temp_vector_3D(2)**2+temp_vector_3D(3)**2)
 
             end do
-            end if
 
             u_taub_squared(i) = u_taub**2
         end do
@@ -1861,7 +1609,7 @@
         if (surface_allocated) then
             do i=1,NNodes_sur
                 temp_vector_1D = node_val(surface_forcing,i)
-                u_taus_squared(i) = max(1e-12,abs(temp_vector_1D(1)))
+                u_taus_squared(i) = max(1e-12,temp_vector_1D(1))
                 !  use the Charnock formula to compute the surface roughness
                 z0s(i)=charnock_val*u_taus_squared(i)/gravity_magnitude
                 if (z0s(i).lt.z0s_min) z0s(i)=z0s_min
@@ -1874,11 +1622,11 @@
 
         do i=1,NNodes_bot
             temp_vector_2D = node_val(bottom_velocity,i)
-            u_taub = sqrt(temp_vector_2D(1)**2+temp_vector_2D(2)**2)
+            u_taub = max(1e-12,sqrt(temp_vector_2D(1)**2+temp_vector_2D(2)**2))
 
             !  iterate bottom roughness length MaxIter times
             do ii=1,MaxIter
-                z0b(i)=(1e-7/(max(1e-6,u_taub)+0.03*0.1))
+                z0b(i)=1e-7/(max(1e-6,u_taub)+0.03*0.1)
                 rr=kappa/log(z0b(i))
 
                 ! compute the friction velocity at the bottom

=== added file 'parameterisation/Ice_Melt_Interface.F90'
--- parameterisation/Ice_Melt_Interface.F90	1970-01-01 00:00:00 +0000
+++ parameterisation/Ice_Melt_Interface.F90	2012-09-04 15:00:37 +0000
@@ -0,0 +1,950 @@
+!    Copyright (C) 2006 Imperial College London and others.
+!
+!    Please see the AUTHORS file in the main source directory for a full list
+!    of copyright holders.
+!
+!    Prof. C Pain
+!    Applied Modelling and Computation Group
+!    Department of Earth Science and Engineering
+!    Imperial College London
+!
+!    [e-mail]
+!
+!    This library is free software; you can redistribute it and/or
+!    modify it under the terms of the GNU Lesser General Public
+!    License as published by the Free Software Foundation,
+!    version 2.1 of the License.
+!
+!    This library is distributed in the hope that it will be useful,
+!    but WITHOUT ANY WARRANTY; without even the implied arranty of
+!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+!    Lesser General Public License for more details.
+!
+!    You should have received a copy of the GNU Lesser General Public
+!    License along with this library; if not, write to the Free Software
+!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+!    USA
+
+#include "fdebug.h"
+
+module ice_melt_interface
+  use quadrature
+  use elements
+  use field_derivatives
+  use fields
+  use fields_base
+  use field_options
+  use fields_allocates
+  use state_module
+  use spud
+  use global_parameters, only: FIELD_NAME_LEN, OPTION_PATH_LEN
+  use state_fields_module
+  use boundary_conditions
+  use fields_manipulation
+  use surface_integrals
+  use fetools
+  use vector_tools
+  use FLDebug
+  use integer_set_module
+  use pickers_inquire
+  use transform_elements
+  use state_fields_module
+  use global_parameters, only: PYTHON_FUNC_LEN
+  use embed_python
+implicit none
+
+  private
+  integer, save                   :: nnodes, dimen
+  type(vector_field), save        :: surface_positions
+  type(vector_field), save        :: funky_positions
+  ! Normals to the iceshelf interface surface
+  type(vector_field), save        :: unit_normal_vectors
+  integer, dimension(:), pointer, save :: sf_nodes_point
+
+  ! Fields and variables for the interface surface
+  type(scalar_field), save        :: ice_surfaceT, ice_surfaceS! these are used to populate the bcs
+
+  public :: melt_interface_initialisation, populate_iceshelf_boundary_conditions, melt_interface_calculate, melt_interface_boundary_condition, melt_interface_cleanup, melt_interface_allocate_surface
+
+contains
+
+  ! Calculation of the (ice) interface meltrate using the parameterisation described in Holland et al. 2008 doi:10.1175/2007JCLI1909.1
+
+  !----------
+  ! initialise parameters based on options
+  !----------
+
+  subroutine melt_interface_initialisation(state)
+
+    type(state_type), intent(inout)     :: state
+    type(scalar_field), pointer         :: T, S
+    ! For the case when Dirichlet boundary conditions are used
+    character(len=FIELD_NAME_LEN)       :: bc_type
+    type(integer_set)                   :: surface_ids
+    integer, dimension(:),allocatable   :: interface_surface_id
+    integer, dimension(2)               :: shape_option
+    type(mesh_type), pointer            :: mesh
+    type(mesh_type)                     :: surface_mesh
+    ! Allocated and returned by create_surface_mesh
+    integer, dimension(:), pointer      :: surface_nodes
+    integer, dimension(:), allocatable  :: surface_element_list
+    integer                             :: i,the_node
+    ! For input of the hydrostatic pressure
+    character(len=PYTHON_FUNC_LEN)      :: subshelf_hydrostatic_pressure_function
+    real                                :: current_time
+    type(vector_field), pointer         :: positions
+    type(scalar_field), pointer         :: subshelf_hydrostatic_pressure
+    real :: c0, cI, L, TI, a, b, gammaT, gammaS, farfield_distance
+    real                                :: T_steady, S_steady
+    logical :: calculate_boundaries_T, calculate_boundaries_S
+    logical :: calculate_boundary_temperature, calculate_boundary_salinity
+    character(len=*), parameter         :: option_path = '/ocean_forcing/iceshelf_meltrate/Holland08'
+
+    ewrite(1,*) 'Melt interface initialisation begins'
+
+    call melt_interface_read_coefficients(c0=c0, cI=cI, L=L, TI=TI, a=a, b=b, gammaT=gammaT, gammaS=gammaS, farfield_distance=farfield_distance, T_steady=T_steady, S_steady=S_steady)
+
+    calculate_boundary_temperature=have_option(trim(option_path)//'/calculate_boundaries/bc_value_temperature')
+    calculate_boundary_salinity=have_option(trim(option_path)//'/calculate_boundaries/bc_value_salinity')
+    if (calculate_boundary_temperature) write(3,*) "Melt interface, initialised T_steady", T_steady
+    if (calculate_boundary_salinity) write(3,*) "Melt interface, initialised S_steady", S_steady
+
+    ! bc= Dirichlet initialize T and S at the ice-ocean interface
+    ! This change with bc type
+    if (have_option(trim(option_path)//'/calculate_boundaries')) then
+      call get_option(trim(option_path)//'/calculate_boundaries', bc_type)
+      select case(bc_type)
+      case('neumann')
+        ewrite(3,*) 'Calculating steady Neumann boundary condition - nothing to be done here.'
+      case('dirichlet')
+        ewrite(3,*) 'Calculating steady Dirichlet boundary condition - nothing to be done here.'
+        ! Define values of temperature and salinity at ice-ocean interface
+        ! Get the surface_id of the ice-ocean interface
+        ! Note this code needs working through - Neumann is definitely the best choice here!
+        shape_option = option_shape(trim(option_path)//"/melt_surfaceID")
+        allocate(interface_surface_id(1:shape_option(1)))
+        call get_option(trim(option_path)//'/melt_surfaceID', interface_surface_id)
+        call allocate(surface_ids)
+        call insert(surface_ids,interface_surface_id)
+        mesh => extract_mesh(state,"VelocityMesh")
+        ! Obtain surface_mesh, surface_element_list, from mesh and surface_id
+        call melt_surf_mesh(mesh,surface_ids,surface_mesh,surface_nodes,surface_element_list)
+
+        T => extract_scalar_field(state,"Temperature")
+        S => extract_scalar_field(state,"Salinity")
+        do i=1, size(surface_nodes)
+          the_node = surface_nodes(i)
+          call set(T,the_node,T_steady)
+          call set(S,the_node,S_steady)
+        enddo
+      case default
+        FLAbort('Unknown boundary type for ice-ocean interface.')
+      end select
+    else
+
+    endif
+
+    ! Read the hydrostatic pressure
+    if(have_option(trim(option_path)//'/subshelf_hydrostatic_pressure')) then
+       call get_option(trim(option_path)//'/subshelf_hydrostatic_pressure/python', subshelf_hydrostatic_pressure_function)
+       ! Get current time
+         call get_option("/timestepping/current_time", current_time)
+       ! Set initial condition from python function
+       ! TODO: This field should be generated automatically and is not needed in the schema
+       subshelf_hydrostatic_pressure => extract_scalar_field(state,"subshelf_hydrostatic_pressure")
+       positions => extract_vector_field(state,"Coordinate")
+       call set_from_python_function(subshelf_hydrostatic_pressure, trim(subshelf_hydrostatic_pressure_function), positions, current_time)
+       ewrite(1,*) "Melt interface initialisation, found hydrostatic pressure field"
+    endif
+
+    ! Additional initialisation routines
+    call melt_interface_allocate_surface(state)
+    call melt_interface_calculate(state)
+    ! Boundary conditions for the (ice) melt interface parameterisation
+    if (have_option(trim(option_path)//'/calculate_boundaries')) then
+      call melt_interface_boundary_condition(state)
+    endif
+
+    ewrite(1,*) "Melt interface initialisation end"
+
+  end subroutine melt_interface_initialisation
+
+  subroutine populate_iceshelf_boundary_conditions(state)
+    type(state_type), intent(in)       :: state
+    type(scalar_field), pointer        :: T,S
+    character(len=FIELD_NAME_LEN)        :: bc_type
+    integer, dimension(:), allocatable   :: surf_id
+    integer, dimension(2)                :: shape_option
+    type(integer_set)                    :: surface_ids
+    !!
+    type(mesh_type), pointer           :: surface_mesh
+    type(scalar_field)                 :: scalar_surface_field
+
+    ewrite(1,*) "-----*** Begin iceshelf BC-----"
+    ! Get vector of surface ids
+    shape_option=option_shape("/ocean_forcing/iceshelf_meltrate/Holland08/melt_surfaceID")
+    allocate(surf_id(1:shape_option(1)))
+    call get_option("/ocean_forcing/iceshelf_meltrate/Holland08/melt_surfaceID", surf_id)
+    ewrite(1,*) "surf_id", surf_id
+    call get_option("/ocean_forcing/iceshelf_meltrate/Holland08/calculate_boundaries", bc_type)
+    call allocate(surface_ids)
+    call insert(surface_ids,surf_id)
+    ewrite(1,*) "set2vector(surface_ids)", set2vector(surface_ids)
+
+    ewrite(1,*) "bc_type: ", bc_type
+     ! Add boundary condition
+    T => extract_scalar_field(state,"Temperature")
+    S => extract_scalar_field(state,"Salinity")
+!    do i=1,node_count(T)
+!        ewrite(1,*) "line2100,i,node_val(T,i): ",i, node_val(T,i)
+!    enddo
+    call add_boundary_condition(T, 'temperature_iceshelf_BC', bc_type, surf_id)
+    call add_boundary_condition(S, 'salinity_iceshelf_BC', bc_type, surf_id)
+    deallocate(surf_id)
+
+    ! mesh of only the part of the surface where this b.c. applies for temperature
+    call get_boundary_condition(T, 'temperature_iceshelf_BC', surface_mesh=surface_mesh)
+    call allocate(scalar_surface_field, surface_mesh, name="value")
+    call insert_surface_field(T, 'temperature_iceshelf_BC', scalar_surface_field)
+    call deallocate(scalar_surface_field)
+
+    ! mesh of only the part of the surface where this b.c. applies for salinity
+    call get_boundary_condition(S, 'salinity_iceshelf_BC', surface_mesh=surface_mesh)
+    call allocate(scalar_surface_field, surface_mesh, name="value")
+    call insert_surface_field(S, 'salinity_iceshelf_BC', scalar_surface_field)
+    call deallocate(scalar_surface_field)
+
+    ewrite(1,*) "-----*** End iceshelf BC-----"
+  end subroutine populate_iceshelf_boundary_conditions
+
+  subroutine melt_interface_calculate(state)
+    ! Calculate the melt rate
+
+    type(state_type), intent(inout)     :: state
+    type(scalar_field), pointer         :: Tb, Sb,MeltRate,Heat_flux,Salt_flux
+    type(scalar_field), pointer         :: scalarfield
+    type(vector_field), pointer         :: velocity,positions
+    ! Debugging purpose pointer
+    type(scalar_field), pointer         :: T_loc,S_loc,P_loc
+    type(vector_field), pointer         :: V_loc,Location,Location_org
+    real, dimension(:), allocatable     :: vel
+    integer                             :: i
+    ! Some internal variables
+    real                                :: speed, T,S,P,Aa,Bb,Cc,topo
+    real :: c0, cI, L, TI, a, b, gammaT, gammaS, farfield_distance
+    real                                ::loc_Tb,loc_Sb,loc_meltrate,loc_heatflux,loc_saltflux
+    ! Aa*Sb^2+Bv*Sb+Cc
+    ! Sink mesh part
+    integer                             :: ele,stat,the_node
+    real, dimension(:), allocatable     :: local_coord,coord
+    integer, dimension(:), allocatable  :: surface_node_list
+    type(scalar_field)                  :: re_pressure,re_DistanceToTop
+    type(scalar_field), pointer         :: temperature, salinity
+
+    ewrite(1,*) "Melt interface calculation begins"
+
+    call melt_interface_read_coefficients(c0=c0, cI=cI, L=L, TI=TI, a=a, b=b, gammaT=gammaT, gammaS=gammaS, farfield_distance=farfield_distance)
+
+    !! All the variable under /ocean_forcing/iceshelf_meltrate/Holland08 should be in coordinate mesh
+    !! coordinate mesh = continous mesh
+    MeltRate => extract_scalar_field(state,"MeltRate")
+    call set(MeltRate,0.0)
+    Tb => extract_scalar_field(state,"Tb")
+    Sb => extract_scalar_field(state,"Sb")
+    ! call set(MeltRate,setnan(arg))
+    call set(Tb,0.0)
+    call set(Sb,-1.0)
+    Heat_flux => extract_scalar_field(state,"Heat_flux")
+    Salt_flux => extract_scalar_field(state,"Salt_flux")
+    call set(Heat_flux,0.0)
+    call set(Salt_flux,0.0)
+
+    ! Local far field values
+    T_loc => extract_scalar_field(state,"Tloc")
+    S_loc => extract_scalar_field(state,"Sloc")
+    P_loc => extract_scalar_field(state,"Ploc")
+    V_loc => extract_vector_field(state,"Vloc")
+    Location => extract_vector_field(state,"Location")
+    Location_org => extract_vector_field(state,"Location_org")
+    call set(T_loc,0.0)
+    call set(S_loc,-1.0)
+    call set(P_loc,-1.0)
+    allocate(vel(V_loc%dim))
+    vel = 0.0
+    call set(V_loc,vel)
+    call set(Location,vel)
+    call set(Location_org,vel)
+
+    positions => extract_vector_field(state,"Coordinate")
+
+    ! Surface node list
+    allocate(surface_node_list(size(sf_nodes_point)))
+    ! sf_nodes is calculated in "melt_allocate_surface"
+    surface_node_list=sf_nodes_point
+
+    ! Pressure read the hydrostatic pressure, specified at schema
+    scalarfield => extract_scalar_field(state,"subshelf_hydrostatic_pressure")
+    call allocate(re_pressure,positions%mesh, name="RePressure")
+    call remap_field(scalarfield,re_pressure,stat)
+
+    if (have_option("/geometry/ocean_boundaries/scalar_field::DistanceToBottom")) then
+        scalarfield => extract_scalar_field(state,"DistanceToBottom")
+        call allocate(re_DistanceToTop,positions%mesh, name="ReDistanceToBottom")
+        call remap_field(scalarfield,re_DistanceToTop,stat)
+    endif
+
+    allocate(local_coord(positions%dim+1))
+    allocate(coord(positions%dim))
+
+    temperature => extract_scalar_field(state,"Temperature")
+    salinity    => extract_scalar_field(state,"Salinity")
+    velocity => extract_vector_field(state,"Velocity")
+
+    ! Loop over the surface nodes to calculate melt rate
+    do i=1,size(surface_node_list)
+        the_node = surface_node_list(i)
+        ! Interpolating
+        coord = node_val(funky_positions,the_node)
+        call picker_inquire(positions, coord, ele, local_coord,global=.false.)
+
+        !! If sum(local_coord) is not equal to 1,
+        !! we know that this coord (funky position) does not exist in the domain.
+        if (sum(local_coord) .gt. 2.0) then
+            ewrite(0,*) "Melt interface, funky coordinate: ", node_val(funky_positions,the_node)
+            !! node_val(surface_positions,the_node) = node_val(positions,the_node)
+            ewrite(0,*) "Melt interface, original element: ",ele
+            ewrite(0,*) "Melt interface, sum of local_coord: ",  sum(local_coord)
+            FLExit("Melt interface, your funky_positions is out of the domain. Change melt_LayerLength.")
+        endif
+
+        !Number of nodes per element for temperature
+        ! Find T at a particular element given the field,element,
+        !INPUT: field, element number
+        !Output: real value
+
+        call scalar_finder_ele(temperature,ele,positions%dim+1,local_coord,T)
+        call scalar_finder_ele(salinity,ele,positions%dim+1,local_coord,S)
+        call vector_finder_ele(velocity,ele,positions%dim+1,local_coord,vel)
+
+        ! Pressure from the schema
+        P = node_val(re_pressure,the_node)
+        speed = sqrt(sum(vel**2))
+
+        if (speed .lt. 0.001) then
+            speed = 0.001
+            !TODO: Use a maxval here instead
+        endif
+        ! constant = -7.53e-8 [C Pa^(-1)] comes from Holland and Jenkins Table 1
+        ! TODO: Define as a constant explicitly, i.e. real, parameter :: topo = -7.53e-8
+        topo = -7.53e-8*P
+
+        ! Define Aa,Bb,Cc
+        ! Aa*Sb**2 + Bb*Sb + Cc = 0.0
+        Aa = -gammaS * speed * cI * a + a * c0 * gammaT * speed
+
+        Bb = -gammaS * speed * L + gammaS * speed * S * cI * a
+        Bb = Bb - gammaS * speed * cI * (b + topo) + gammaS * speed * cI * TI
+        Bb = Bb - c0 * gammaT * speed * T + c0 * gammaT * speed * (b + topo)
+
+        Cc = gammaS * speed * S * L + gammaS * speed * S * cI * (b + topo) + gammaS * speed * S * (-cI * TI)
+
+        ! This could be a linear equation if Aa=0
+        if (Aa .eq. 0.0) then
+            loc_Sb = -Cc/Bb
+        else
+          ! Calculate for the 2nd oewrite(1,*) "size(surface_element_list)"rder polynomial.
+          ! We have two solutions.
+          loc_Sb = (-Bb + sqrt(Bb**2 - 4.0*Aa*Cc))/(2.0*Aa)
+          ! loc_Sb has to be larger than 0; since the salinity in the ocean is positive definite.
+          if (loc_Sb .lt. 0.0) then
+            loc_Sb = (-Bb - sqrt(Bb**2 - 4.0*Aa*Cc))/(2.0*Aa)
+          endif
+          if (loc_Sb .lt. 0.0) then
+            ewrite(0,*) "Melt interface, loc_Sb: ",  loc_Sb
+            ewrite(0,*) "Melt interface, Aa: ",  Aa
+            ewrite(0,*) "Melt interface, Bb: ",  Bb
+            ewrite(0,*) "Melt interface, Cc: ",  Cc
+            ewrite(0,*) "Melt interface, T: ",  T
+            ewrite(0,*) "Melt interface, S: ",  S
+            ewrite(0,*) "Melt interface, P: ",  P
+            ewrite(0,*) "Melt interface, speed: ",  speed
+            FLExit("Melt interface, Sb is negative. The range of Salinity is not right.")
+          endif
+        endif
+
+        loc_Tb = a*loc_Sb + b + topo
+        loc_meltrate = gammaS*speed*(S-loc_Sb)/loc_Sb
+        !! Heat flux to the ocean
+        loc_heatflux = (gammaT*speed+ loc_meltrate)*(loc_Tb - T) ! or loc_meltrate*L + loc_meltrate*cI*(loc_Tb-TI)
+        !! Salt flux to the ocean
+        loc_saltflux = (gammaS*speed+loc_meltrate)*(loc_Sb - S)
+
+        !! These are needed to implement BCs.
+        call set(MeltRate, the_node, loc_meltrate)
+        call set(Tb, the_node, loc_Tb)
+        call set(Sb, the_node, loc_Sb)
+        call set(Heat_flux, the_node,loc_heatflux)
+        call set(Salt_flux, the_node,loc_saltflux)
+
+        !!Far field values
+        call set(T_loc,the_node,T)
+        call set(S_loc,the_node,S)
+        call set(P_loc,the_node,P)
+        call set(V_loc,the_node,vel)
+        call set(Location,the_node,node_val(funky_positions,the_node))
+        call set(Location_org,the_node,node_val(positions,the_node))
+
+    enddo
+
+    deallocate(local_coord)
+    deallocate(coord)
+    call deallocate(re_pressure)
+    ! TODO: check this below
+    if (have_option("/geometry/ocean_boundaries/scalar_field::DistanceToBottom")) then
+        call deallocate(re_DistanceToTop)
+    endif
+
+    ! Remap meltrate_p heat_flux_p salt_flux_p onto MeltRate Heat_flux Salt_flux
+    ! Mappting continous mesh to discontinous mesh
+
+    ewrite(1,*) "Melt interface calculation end"
+
+  end subroutine melt_interface_calculate
+
+  subroutine melt_interface_boundary_condition(state)
+    ! See preprocessor/Boundary_Conditions_From_options.F90 populate_iceshelf_boundary_conditions(states(1)) as well
+    type(state_type), intent(inout)     :: state
+    ! Boundary conditions
+    type(scalar_field), pointer         :: TT,SS
+    type(scalar_field), pointer         :: scalar_surface, scalar_surface2
+    type(mesh_type), pointer            :: ice_mesh
+    character(len=FIELD_NAME_LEN)       :: bc_type
+    type(scalar_field)                  :: T_bc,S_bc
+    type(scalar_field), pointer         :: Tflux,Sflux
+    type(scalar_field), pointer         :: Tb,Sb
+    integer                             :: i, the_node
+    real                                :: Tz,Sz
+    type(mesh_type), pointer            :: mesh
+    type(mesh_type)                     :: surface_mesh
+    integer, dimension(:), pointer      :: surface_nodes ! allocated and returned by create_surface_mesh
+    integer, dimension(:), allocatable  :: surface_element_list
+    type(integer_set)                   :: surface_ids
+    integer, dimension(:),allocatable   :: surf_id
+    integer, dimension(2) :: shape_option
+    type(vector_field)                   :: unit_normal_vectors_vel
+    type(scalar_field)                  :: heat_flux_vel,salt_flux_vel
+    type(vector_field), pointer         :: velocity
+    integer                             :: stat
+    real, dimension(:), allocatable     :: vel
+    real                                :: T_steady, S_steady
+    logical :: calculate_boundary_temperature, calculate_boundary_salinity
+    character(len=*), parameter         :: option_path = '/ocean_forcing/iceshelf_meltrate/Holland08'
+    character(len=*), parameter         :: option_path_bc = '/ocean_forcing/iceshelf_meltrate/Holland08/calculate_boundaries'
+
+    ewrite(1,*) "Melt interface boundary condition begins"
+    call melt_interface_read_coefficients(T_steady=T_steady, S_steady=S_steady)
+    calculate_boundary_temperature=have_option(trim(option_path)//'/calculate_boundaries/bc_value_temperature')
+    calculate_boundary_salinity=have_option(trim(option_path)//'/calculate_boundaries/bc_value_salinity')
+
+!! See ./preprocessor/Boundary_Conditions_From_options.F90 populate_iceshelf_boundary_conditions(states(1)) as well
+    ! Get the surface_id of the ice-ocean interface
+    shape_option=option_shape(trim(option_path)//"/melt_surfaceID")
+    allocate(surf_id(1:shape_option(1)))
+    call get_option(trim(option_path)//'/melt_surfaceID',surf_id)
+    call allocate(surface_ids)
+    call insert(surface_ids,surf_id)
+
+    TT=> extract_scalar_field(state,"Temperature")
+    SS=> extract_scalar_field(state,"Salinity")
+    velocity => extract_vector_field(state,"Velocity")
+
+    if (node_count(TT) .eq. node_count(velocity)) then
+      mesh => extract_mesh(state,"VelocityMesh")
+    else
+      mesh => extract_mesh(state,"CoordinateMesh")
+    endif
+    call melt_surf_mesh(mesh,surface_ids,surface_mesh,surface_nodes,surface_element_list)
+
+!! Remapt the contious mesh to discon, unit_normal vecotr
+
+    call allocate(unit_normal_vectors_vel,velocity%dim,TT%mesh,"MyVelocitySurfaceMesh")
+
+    allocate(vel(velocity%dim))
+    vel = 0.0
+    call set(unit_normal_vectors_vel,vel)
+    deallocate(vel)
+
+    call remap_field(unit_normal_vectors,unit_normal_vectors_vel,stat)
+
+    call allocate(T_bc,TT%mesh, name="T_boundary")
+    call allocate(S_bc,SS%mesh, name="S_boundary")
+
+    ! Surface node list
+    ! This change with bc type
+    call get_option(trim(option_path_bc), bc_type)
+
+    select case(bc_type)
+    case("neumann")
+      Tflux => extract_scalar_field(state,"Heat_flux")
+      Sflux => extract_scalar_field(state,"Salt_flux")
+
+      !! Remap the heat flux
+      call allocate(heat_flux_vel,TT%mesh,"MyHeatFluxSurfaceMesh")
+      call set(heat_flux_vel,0.0)
+      call remap_field(Tflux,heat_flux_vel,stat)
+
+      !! Remap the salt flux
+      call allocate(salt_flux_vel,SS%mesh,"MySaltFluxSurfaceMesh")
+      call set(salt_flux_vel,0.0)
+      call remap_field(Sflux,salt_flux_vel,stat)
+
+      ! create a surface mesh to place values onto. This is for the top surface
+      call get_boundary_condition(TT, 'temperature_iceshelf_BC', surface_mesh=ice_mesh)
+      call allocate(ice_surfaceT, ice_mesh, name="ice_surfaceT")
+      !! Temperature
+      scalar_surface => extract_surface_field(TT, 'temperature_iceshelf_BC', "value")
+
+      do i=1,node_count(ice_surfaceT)
+        the_node = surface_nodes(i)
+        !                Kt = (sum(node_val(unit_normal_vectors_vel,the_node)*Kt_ar))
+        !                Tz = node_val(heat_flux_vel,the_node)/(-Kt*c0)
+        Tz = node_val(heat_flux_vel,the_node)
+        call set(scalar_surface,i,Tz)
+      enddo
+
+      ! Salinity
+      call get_boundary_condition(SS, 'salinity_iceshelf_BC', surface_mesh=ice_mesh)
+      call allocate(ice_surfaceS, ice_mesh, name="ice_surfaceS")
+      !! Salinity
+      scalar_surface2 => extract_surface_field(SS, 'salinity_iceshelf_BC', "value")
+      do i=1,node_count(ice_surfaceS)
+        the_node = surface_nodes(i)
+        !                Ks = abs(sum(node_val(unit_normal_vectors_vel,the_node)*Ks_ar))
+        !                Sz = node_val(salt_flux_vel,the_node)/(-Ks)
+        Sz = node_val(salt_flux_vel,the_node)
+        !call set(ice_surfaceS,i,Sz)
+        call set(scalar_surface2,i,Sz)
+      enddo
+      ! Deallocate the heat_flux and salt_flux on velocity mesh
+      call deallocate(heat_flux_vel)
+      call deallocate(salt_flux_vel)
+
+      ! If the fluxes were constant
+      if (calculate_boundary_temperature) then
+        call set(T_bc,T_steady)
+        ! create a surface mesh to place values onto. This is for the top surface
+        call get_boundary_condition(TT, 'temperature_iceshelf_BC', surface_mesh=ice_mesh)
+        call allocate(ice_surfaceT, ice_mesh, name="ice_surfaceT")
+        do i=1,node_count(ice_surfaceT)
+          the_node = surface_nodes(i)
+          call set(ice_surfaceT,i,node_val(T_bc,the_node))
+        enddo
+        ! Temperature
+        scalar_surface => extract_surface_field(TT, 'temperature_iceshelf_BC', "value")
+        call remap_field(ice_surfaceT, scalar_surface)
+      endif
+
+      if (have_option(trim(option_path_bc)//'/bc_value_salinity')) then
+        call set(S_bc,S_steady)
+        ! Salinity
+        call get_boundary_condition(SS, 'salinity_iceshelf_BC', surface_mesh=ice_mesh)
+        call allocate(ice_surfaceS, ice_mesh, name="ice_surfaceS")
+        do i=1,node_count(ice_surfaceS)
+          the_node = surface_nodes(i)
+          call set(ice_surfaceS,i,node_val(S_bc,the_node))
+        enddo
+        !! Salinity
+          scalar_surface => extract_surface_field(SS, 'salinity_iceshelf_BC', "value")
+          call remap_field(ice_surfaceS, scalar_surface)
+      endif
+
+    case("dirichlet")
+      Tb => extract_scalar_field(state,"Tb")
+      Sb => extract_scalar_field(state,"Sb")
+      do i=1,node_count(T_bc)
+        call set(T_bc,i,node_val(Tb,i))
+        call set(S_bc,i,node_val(Sb,i))
+      enddo
+      ! When testing BC implementation
+      if (calculate_boundary_temperature) then
+        call set(T_bc,T_steady)
+        ewrite(3,*) "Melt interface, initialised T_steady", T_steady
+      endif
+      if (calculate_boundary_salinity) then
+        call set(S_bc,S_steady)
+        ewrite(3,*) "Melt interface, initialised S_steady", S_steady
+      endif
+
+    case default
+      FLAbort('Unknown boundary condition for ice-ocean interface')
+    end select
+
+    ! create a surface mesh to place values onto. This is for the top surface
+    !call get_boundary_condition(TT, 'temperature_iceshelf_BC', surface_mesh=ice_mesh)
+    !call allocate(ice_surfaceT, ice_mesh, name="ice_surfaceT")
+    ! Define ice_surfaceT according to Heat_flux?
+
+    !do i=1,node_count(ice_surfaceT)
+    !    the_node = surface_node_list(i)
+    !    call set(ice_surfaceT,i,node_val(T_bc,the_node))
+    !enddo
+
+    ! Salinity
+    !call get_boundary_condition(SS, 'salinity_iceshelf_BC', surface_mesh=ice_mesh)
+    !call allocate(ice_surfaceS, ice_mesh, name="ice_surfaceS")
+    !do i=1,node_count(ice_surfaceS)
+    !    the_node = surface_node_list(i)
+    !    call set(ice_surfaceS,i,node_val(S_bc,the_node))
+    !enddo
+
+    !!! Temperature
+    !scalar_surface => extract_surface_field(TT, 'temperature_iceshelf_BC', "value")
+    !call remap_field(ice_surfaceT, scalar_surface)
+    !!! Salinity
+    !scalar_surface => extract_surface_field(SS, 'salinity_iceshelf_BC', "value")
+    !call remap_field(ice_surfaceS, scalar_surface)
+
+    call deallocate(T_bc)
+    call deallocate(S_bc)
+    call deallocate(ice_surfaceT)
+    call deallocate(ice_surfaceS)
+    call deallocate(unit_normal_vectors_vel)
+    call deallocate(surface_mesh)
+
+    ewrite(1,*) "Melt interface boundary condition end"
+
+  end subroutine melt_interface_boundary_condition
+
+
+  subroutine melt_interface_cleanup()
+    ewrite(1,*) "Melt interface, cleaning up melt variables"
+    call deallocate(surface_positions)
+    call deallocate(funky_positions)
+    call deallocate(unit_normal_vectors)
+
+    deallocate(sf_nodes_point)
+  end subroutine melt_interface_cleanup
+
+
+  !------------------------------------------------------------------!
+  !                       Private subroutines                        !
+  !------------------------------------------------------------------!
+
+  subroutine melt_interface_allocate_surface(state)
+    type(state_type), intent(inout)     :: state
+    type(vector_field), pointer         :: positions
+    type(mesh_type), pointer            :: mesh
+    type(mesh_type)                     :: surface_mesh
+    ! Allocated and returned by create_surface_mesh
+    integer, dimension(:), pointer      :: surface_nodes
+    integer, dimension(:), allocatable  :: surface_element_list
+    integer    :: face,i,j,k
+    real, dimension(:), allocatable     :: coord
+    ! For transform_facet_to_physical
+    real, dimension(:,:), allocatable   :: normal
+    real, dimension(:), allocatable     :: av_normal !Average of normal
+    type(element_type), pointer     :: x_shape_f
+    real, dimension(:), allocatable     :: xyz !New location of surface_mesh, farfield_distance away from the boundary
+    type(integer_set)                   :: surface_ids
+    integer, dimension(:),allocatable   :: interface_surface_id
+
+    ! For the vector averaging schem, taking the area of the surface element etc
+    real, dimension(:,:), allocatable    :: table
+    integer, dimension(:,:), allocatable :: node_occupants
+    integer                              :: st,en,node,dim_vec
+    real                                 :: area_sum
+    integer, dimension(2) :: shape_option
+    real, dimension(:), allocatable     :: local_coord
+    integer                             :: ele, counter,Nit
+    real, dimension(:), allocatable     :: vel
+    real :: c0, cI, L, TI, a, b, gammaT, gammaS, farfield_distance
+    character(len=*), parameter         :: option_path = '/ocean_forcing/iceshelf_meltrate/Holland08/melt_surfaceID'
+
+    ewrite(1,*) "Melt interface allocation of surface parameters"
+
+    call melt_interface_read_coefficients(c0=c0, cI=cI, L=L, TI=TI, a=a, b=b, gammaT=gammaT, gammaS=gammaS, farfield_distance=farfield_distance)
+
+    ! Get the surface_id of the ice-ocean interface
+    shape_option=option_shape(trim(option_path))
+    allocate(interface_surface_id(1:shape_option(1)))
+    call get_option(trim(option_path), interface_surface_id)
+    call allocate(surface_ids)
+    call insert(surface_ids,interface_surface_id)
+
+    mesh => extract_mesh(state,"CoordinateMesh")
+    ! Input, mesh, surface_id
+    ! Output surface_mesh, surface_element_list
+    call melt_surf_mesh(mesh, surface_ids, surface_mesh, surface_nodes, surface_element_list)
+
+    allocate(sf_nodes_point(size(surface_nodes)))
+    sf_nodes_point = surface_nodes
+    positions => extract_vector_field(state,"Coordinate")
+    ! Sink the surface_mesh to the 'far field'
+    call allocate(surface_positions,positions%dim,mesh,"MySurfacePosition")
+    call allocate(funky_positions,surface_positions%dim, surface_positions%mesh,"MyFunkyPosition")
+    call allocate(unit_normal_vectors,surface_positions%dim, surface_positions%mesh,"MyFunkyUnitVec")
+    allocate(vel(positions%dim))
+    vel = 0.0
+    call set(unit_normal_vectors,vel)
+    deallocate(vel)
+    ! Check above
+
+    ! Remap the positions vector to the surface mesh, surface_positions
+    call remap_field_to_surface(positions, surface_positions, surface_element_list)
+
+    ! Loop over # of surface elements
+    !nf = size(surface_element_list)
+    !2 comes because there are 2 nodes per a surface_element (assumption!)
+    ! this number could be three for 3D problem, room for change
+    ! dim_vec is the dimension of the positions, either 2 or 3.
+    ! node_occupants(1,1:dim_vec*nf) = nodes that occupies the surface elements
+    ! node_occupants(2,:) = surface elements
+    ! table(1,:) = area of the face
+    ! table(2,:) = normal_x
+    ! table(3,:) = normal_y
+    ! table(4,:) = normal_z
+    dim_vec = positions%dim
+    allocate(coord(dim_vec))
+    allocate(node_occupants(2,dim_vec*size(surface_element_list)))
+    allocate(table(dim_vec+1,dim_vec*size(surface_element_list)))
+    allocate(av_normal(dim_vec))
+    allocate(xyz(dim_vec))
+    allocate(local_coord(positions%dim+1))
+
+    do i=1,size(surface_element_list)
+        st = 1+dim_vec*(i-1)
+        en = dim_vec+dim_vec*(i-1)
+        face = surface_element_list(i)
+        node_occupants(2,st:en) = face
+        node_occupants(1,st:en) = face_global_nodes(positions, face)
+        ! Calculate the area of the surface element
+        ! For 2D, the area is the length
+        if (dim_vec .eq. 2) then
+            table(1,st:en) = calc_area2(positions,node_occupants(1,st:en))
+        endif
+
+        ! For 3D, the area become the area of the triangle (assumption here).
+        ! Subroutine calc_area3 uses Heron's formula.
+        if (dim_vec .eq. 3) then
+            table(1,st:en) = calc_area3(positions,node_occupants(1,st:en))
+        endif
+        ! Compute the normal vector at the surface element.
+        x_shape_f=>face_shape(positions,face)
+        allocate(normal(dim_vec,x_shape_f%ngi)) !(dim x x_shape_f%ngi)
+        call transform_facet_to_physical(positions, face, normal=normal)
+        av_normal = sum(normal,2)
+        !"normal" should be the same; since the normal vector on the quadrature point does not change.
+        av_normal = av_normal/(sum(av_normal*av_normal))**(0.5) ! normalize the vector
+        deallocate(normal)
+        ! Since av_normal is outward to the boundary, we will put -sign. We want the vector into the boundary
+        av_normal = -av_normal
+        ! Store av_normal, the normal vector averaged over quadrature points, in table
+        do j=1,dim_vec
+            table(j+1,st:en)=av_normal(j)
+        enddo
+    enddo
+    ! Now loop over surface_nodes
+    !        ewrite(1,*) "table(1,1:3): ", table(1,1:3)
+    !        ewrite(1,*) "table(2,1:3),normal_x: ", table(2,1:3)
+    !        ewrite(1,*) "table(3,1:3),normal_y: ", table(3,1:3)
+    !        ewrite(1,*) "table(4,1:3),normal_z: ", table(4,1:3)
+    ! In this loop, we will average adjacent normal vectors, using the area of the surface elements.
+
+
+    do i=1,size(node_occupants(1,:))
+        node = node_occupants(1,i)
+        av_normal(:) = 0.0
+        area_sum = 0.0
+        ! This loop average the normal vector based on the areas of surface elements
+        ! which share the "node". Using hash table may speed up this process
+        do j=1,size(node_occupants(1,:))
+            !Pick the surface elements that sheare the "node"
+            if (node_occupants(1,j) .eq. node) then
+                do k=1,dim_vec
+                    !table(1,j) = area of the surface element
+                    av_normal(k) = av_normal(k) + table(k+1,j)*table(1,j)
+                enddo
+                !The total areas of the surface elements that occupies the "node"
+                area_sum = area_sum + table(1,j)
+            endif
+        enddo
+
+        av_normal = av_normal / area_sum
+        ! normalize
+        av_normal = av_normal/(sum(av_normal*av_normal))**(0.5)
+
+        farfield_distance = abs(farfield_distance)
+        ! Shift the location of the surface nodes.
+        !The coordinate of the surface node.
+
+        coord = node_val(positions,node)
+        ! farfield_distance = ||xyz - coord|| xyz and coord are vectors
+        xyz = av_normal*farfield_distance + coord
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+        ! have options /ocean_forcing/iceshelf_meltrate/Holland08/melt_LayerRelax
+        if (have_option("/ocean_forcing/iceshelf_meltrate/Holland08/melt_LayerRelax")) then
+          call get_option('/ocean_forcing/iceshelf_meltrate/Holland08/melt_LayerRelax/N',Nit)
+          !!Check if xyz is within the domain
+          call picker_inquire(positions, xyz, ele, local_coord,global=.false.)
+          !! If sum(local_coord) is not equal to 1,
+          !! we know that this coord (funky position) does not exist in the domain.
+
+          counter = 1
+          do while (sum(local_coord) .gt. 2.0 .and. counter .lt. Nit)
+            xyz = av_normal*(farfield_distance-counter*(farfield_distance/Nit)) + coord
+            call picker_inquire(positions, xyz, ele, local_coord,global=.false.)
+            counter = counter+1
+          enddo
+
+          if (sum(local_coord) .gt. 2.0) then
+            ewrite(1,*) "icehslef location, coord: ", coord
+            ewrite(1,*) "funky position, xyz: ", xyz
+            call picker_inquire(positions, coord, ele, local_coord,global=.false.)
+            !! node_val(surface_positions,the_node) = node_val(positions,the_node)
+            ewrite(1,*) "Original ele: ",ele
+            ewrite(1,*) "sum of local_coord: ",  sum(local_coord)
+            FLExit("In setting funky_position, your funky_positions is out of the domain. Change melt_LayerLength.")
+          endif
+        endif
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+        call set(funky_positions,node,xyz) ! Set the coordinate of sinked nodes, funky positions.
+        call set(surface_positions,node,coord) !Original coordinate of the surface
+        call set(unit_normal_vectors,node,av_normal) !save the unit normal vector for Kt and Ks calculation
+
+        node = 0
+    enddo
+
+    ! Now deallocate junk
+    deallocate(coord)
+    deallocate(table)
+    deallocate(node_occupants)
+    deallocate(av_normal)
+    deallocate(xyz)
+    deallocate(local_coord)
+    call deallocate(surface_ids)
+    call deallocate(surface_mesh)
+    ewrite(1,*) "Melt interface allocation of surface parameters end"
+
+  end subroutine melt_interface_allocate_surface
+
+  subroutine melt_surf_mesh(mesh,surface_ids,surface_mesh,surface_nodes,surface_element_list)
+
+    type(mesh_type), intent(in)                       :: mesh
+    type(integer_set), intent(in)                     :: surface_ids
+    type(mesh_type), intent(out)                      :: surface_mesh
+    integer, dimension(:), pointer, intent(out)       :: surface_nodes ! allocated and returned by create_surface_mesh
+    integer, dimension(:), allocatable, intent(out)   :: surface_element_list
+    !    integer, dimension(:), allocatable                :: surface_nodes_out
+    type(integer_set)                                 :: surface_elements
+    integer :: i
+    ! create a set of surface elements that have surface id in 'surface_ids'
+
+    call allocate(surface_elements)
+    do i=1, surface_element_count(mesh)
+      ! ewrite(1,*) "surf_normal surface_element_id(mesh, i)", surface_element_id(mesh, i)
+      ! ewrite(1,*) "surf_normal surface_ids", set2vector(surface_ids)
+      if (has_value(surface_ids, surface_element_id(mesh, i))) then
+        call insert(surface_elements, i)
+      end if
+    end do
+
+    allocate(surface_element_list(key_count(surface_elements)))
+    surface_element_list=set2vector(surface_elements)
+    call create_surface_mesh(surface_mesh, surface_nodes, mesh, surface_element_list, name=trim(mesh%name)//"ToshisMesh")
+  end subroutine melt_surf_mesh
+
+  subroutine melt_interface_read_coefficients(c0, cI, L, TI, a, b, gammaT, gammaS, farfield_distance, T_steady, S_steady)
+    real, intent(out), optional :: c0, cI, L, TI, a, b, gammaT, gammaS, farfield_distance, T_steady, S_steady
+    real :: Cd
+    character(len=*), parameter :: option_path = '/ocean_forcing/iceshelf_meltrate/Holland08'
+
+    ! Get the 6 model constants
+    ! TODO: Check these exist first and error with a useful message if not - in preprocessor
+    if (present(c0)) call get_option(trim(option_path)//'c0', c0, default = 3974.0)
+    if (present(cI)) call get_option(trim(option_path)//'cI', cI, default = 2009.0)
+    if (present(L))  call get_option(trim(option_path)//'/L', L, default = 3.35e5)
+    if (present(TI)) call get_option(trim(option_path)//'/TI', TI, default = -25.0)
+    if (present(a))  call get_option(trim(option_path)//'/a', a, default = -0.0573)
+    if (present(b))  call get_option(trim(option_path)//'/b', b, default = 0.0832)
+    if (present(farfield_distance)) call get_option(trim(option_path)//'/melt_LayerLength', farfield_distance)
+
+    if (present(gammaT).or.present(gammaS)) call get_option(trim(option_path)//'/Cd', Cd, default = 1.5e-3)
+    if (present(gammaT)) gammaT = sqrt(Cd)/(12.5*(7.0**(2.0/3.0))-9.0)
+    if (present(gammaS)) gammaS = sqrt(Cd)/(12.5*(700.0**(2.0/3.0))-9.0)
+
+    ! When steady boundary options are enabled
+    if (present(T_steady)) call get_option(trim(option_path)//'/calculate_boundaries/bc_value_temperature',T_steady, default = 0.0)
+    if (present(S_steady)) call get_option(trim(option_path)//'/calculate_boundaries/bc_value_salinity',S_steady, default = 34.0)
+
+  end subroutine melt_interface_read_coefficients
+
+  subroutine scalar_finder_ele(scalar,ele,element_dim,local_coord,scalar_out)
+    type(scalar_field), pointer, intent(in)           :: scalar
+    integer, intent(in)                      :: ele
+    integer, intent(in)                       :: element_dim
+    real, dimension(element_dim), intent(in) :: local_coord
+    real, intent(out)                        :: scalar_out
+    integer                        :: j,node
+    integer, dimension(:), allocatable  :: node_lists
+
+    allocate(node_lists(size(ele_nodes(scalar,ele))))
+    node_lists =  ele_nodes(scalar,ele) !Lists of nodes that occupies the element.
+    scalar_out = 0.0
+    do j=1,size(node_lists)
+        node = node_lists(j)
+        scalar_out = scalar_out + node_val(scalar,node)*local_coord(j)
+    enddo
+  end subroutine scalar_finder_ele
+
+  subroutine vector_finder_ele(vector,ele,element_dim,local_coord,vector_out)
+    type(vector_field), pointer, intent(in)           :: vector
+    integer, intent(in)                      :: ele
+    integer, intent(in)                      :: element_dim
+    real, dimension(element_dim), intent(in) :: local_coord
+    real, dimension(element_dim-1), intent(out) :: vector_out
+    integer                        :: j,node
+    integer, dimension(:), allocatable  :: node_lists
+
+    allocate(node_lists(size(ele_nodes(vector,ele))))
+
+    node_lists =  ele_nodes(vector,ele) !Lists of nodes that occupies the element.
+    vector_out(:) = 0.0
+    do j=1,size(node_lists)
+        node = node_lists(j)
+        vector_out = vector_out + node_val(vector,node)*local_coord(j)
+    enddo
+  end subroutine vector_finder_ele
+
+  real function setnan(arg)
+    real :: arg
+    setnan = sqrt(arg)
+  end function
+
+  real function calc_area2(positions,nodes)
+    type(vector_field), pointer        :: positions
+    integer, dimension(2)              :: nodes
+    real, dimension(2)                 :: coord1,coord2
+    coord1 = node_val(positions,nodes(1))
+    coord2 = node_val(positions,nodes(2))
+    calc_area2 = (sum((coord2 - coord1)**2))**0.5
+  end function
+
+  real function calc_area3(positions,nodes)
+    type(vector_field), pointer        :: positions
+    integer, dimension(3)              :: nodes
+    real, dimension(3)                 :: coord1,coord2,coord3
+    real                               :: a,b,c,s
+    coord1 = node_val(positions,nodes(1))
+    coord2 = node_val(positions,nodes(2))
+    coord3 = node_val(positions,nodes(3))
+    ! Use Heron's formula to calculate the area of triangle
+    a = (sum((coord2 - coord1)**2))**0.5
+    b = (sum((coord3 - coord1)**2))**0.5
+    c = (sum((coord2 - coord3)**2))**0.5
+    s = 0.5*(a+b+c)
+    calc_area3 = (s*(s-a)*(s-b)*(s-c))**0.5
+  end function
+
+end module ice_melt_interface

=== modified file 'parameterisation/Makefile.dependencies'
--- parameterisation/Makefile.dependencies	2012-08-29 10:53:27 +0000
+++ parameterisation/Makefile.dependencies	2012-09-04 15:00:37 +0000
@@ -18,16 +18,15 @@
    ../include/fields.mod ../include/fldebug.mod \
    ../include/global_parameters.mod ../include/quadrature.mod \
    ../include/sparse_matrices_fields.mod ../include/spud.mod \
-   ../include/state_fields_module.mod ../include/state_module.mod \
-   ../include/vertical_extrapolation_module.mod
+   ../include/state_fields_module.mod ../include/state_module.mod
 
-../include/iceshelf_meltrate_surf_normal.mod: iceshelf_meltrate_surf_normal.o
+../include/ice_melt_interface.mod: Ice_Melt_Interface.o
 	@true
 
-iceshelf_meltrate_surf_normal.o ../include/iceshelf_meltrate_surf_normal.mod: \
-   iceshelf_meltrate_surf_normal.F90 ../include/boundary_conditions.mod \
+Ice_Melt_Interface.o ../include/ice_melt_interface.mod: \
+   Ice_Melt_Interface.F90 ../include/boundary_conditions.mod \
    ../include/elements.mod ../include/fdebug.h ../include/fetools.mod \
-   ../include/field_derivatives.mod ../include/field_options.mod \
+   ../include/field_derivatives.mod ../include/fields_base.mod ../include/field_options.mod \
    ../include/fields.mod ../include/fields_allocates.mod \
    ../include/fields_manipulation.mod ../include/fldebug.mod \
    ../include/global_parameters.mod ../include/integer_set_module.mod \

=== modified file 'parameterisation/Makefile.in'
--- parameterisation/Makefile.in	2011-10-21 10:05:25 +0000
+++ parameterisation/Makefile.in	2012-09-04 15:00:37 +0000
@@ -47,7 +47,7 @@
 OBJS = Equation_of_State.o \
 Gls_vertical_turbulence_model.o \
 k_epsilon.o \
-iceshelf_meltrate_surf_normal.o
+Ice_Melt_Interface.o
 
 .SUFFIXES: .F90 .c .o .a
 

=== removed file 'parameterisation/iceshelf_meltrate_surf_normal.F90'
--- parameterisation/iceshelf_meltrate_surf_normal.F90	2011-04-07 12:42:58 +0000
+++ parameterisation/iceshelf_meltrate_surf_normal.F90	1970-01-01 00:00:00 +0000
@@ -1,678 +0,0 @@
-!    Copyright (C) 2006 Imperial College London and others.
-!    
-!    Please see the AUTHORS file in the main source directory for a full list
-!    of copyright holders.
-!
-!    Prof. C Pain
-!    Applied Modelling and Computation Group
-!    Department of Earth Science and Engineering
-!    Imperial College London
-!
-!    [e-mail]
-!    
-!    This library is free software; you can redistribute it and/or
-!    modify it under the terms of the GNU Lesser General Public
-!    License as published by the Free Software Foundation,
-!    version 2.1 of the License.
-!
-!    This library is distributed in the hope that it will be useful,
-!    but WITHOUT ANY WARRANTY; without even the implied arranty of
-!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-!    Lesser General Public License for more details.
-!
-!    You should have received a copy of the GNU Lesser General Public
-!    License along with this library; if not, write to the Free Software
-!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
-!    USA
-
-#include "fdebug.h"
-
-module iceshelf_meltrate_surf_normal
-  use quadrature
-  use elements
-  use field_derivatives
-  use fields
-  use field_options
-  use fields_allocates
-  use state_module
-  use spud
-  use global_parameters, only: FIELD_NAME_LEN, OPTION_PATH_LEN
-  use state_fields_module
-  use boundary_conditions
-  use fields_manipulation
-  use surface_integrals
-  use fetools
-  use vector_tools
-  use FLDebug
-  use integer_set_module
-  use pickers_inquire
-  use transform_elements
-  use state_fields_module
-implicit none
-
-  private
-  real, save                      :: c0, cI, L, TI, &
-                                     a, b, gammaT, gammaS,Cd, dist_meltrate
-  integer, save                   :: nnodes,dimen
- 
-  type(vector_field), save        :: surface_positions
-  type(vector_field), save        :: funky_positions 
-  type(integer_set), save         :: sf_nodes !Nodes at the surface
-  !BC
-  integer, dimension(:), pointer, save :: ice_element_list
-  ! these are the fields and variables for the surface values
-  type(scalar_field), save             :: ice_surfaceT,ice_surfaceS! these are used to populate the bcs
-  public :: melt_surf_init, melt_allocate_surface, melt_surf_calc, melt_bc
-
-
-
-contains
-
-!----------
-! initialise parameters based on options
-!----------
-
- subroutine melt_surf_init(state)
-
-    type(state_type), intent(inout)     :: state
-    character(len=OPTION_PATH_LEN)      :: melt_path
-    ! hack
-    type(scalar_field), pointer                  :: T,S
-    ! When bc=Dirichlet 
-     character(len=FIELD_NAME_LEN)       :: bc_type
-    type(integer_set)                   :: surface_ids
-    integer, dimension(:),allocatable   :: surf_id
-    integer, dimension(:), allocatable  :: sf_nodes_ar
-    integer, dimension(2) :: shape_option
-    type(mesh_type), pointer            :: mesh
-    type(mesh_type)                     :: surface_mesh
-    integer, dimension(:), pointer      :: surface_nodes ! allocated and returned by create_surface_mesh
-    integer, dimension(:), allocatable  :: surface_element_list
-     integer                             :: i,the_node
-      
-    melt_path = "/ocean_forcing/iceshelf_meltrate/Holland08"
-
- 
-    ewrite(1,*) "--------Begin melt_init-------------"
-
-    ! Get the 6 model constants
-    call get_option(trim(melt_path)//'/c0', c0, default = 3974.0)
-    call get_option(trim(melt_path)//'/cI', cI, default = 2009.0)
-    call get_option(trim(melt_path)//'/L', L, default = 3.35e5)
-    call get_option(trim(melt_path)//'/TI', TI, default = -25.0)
-    call get_option(trim(melt_path)//'/a', a, default = -0.0573)
-    call get_option(trim(melt_path)//'/b', b, default = 0.0832)
-    call get_option(trim(melt_path)//'/Cd', Cd, default = 1.5e-3)
-    call get_option(trim(melt_path)//'/melt_LayerLength', dist_meltrate)
- 
-    gammaT = sqrt(Cd)/(12.5*(7.0**(2.0/3.0))-9.0)
-    gammaS = sqrt(Cd)/(12.5*(700.0**(2.0/3.0))-9.0)
-
-    !! bc= Dirichlet initialize T and S at the ice-ocean interface
-    !hack
-    !This change with bc type
-    call get_option("/ocean_forcing/iceshelf_meltrate/Holland08/calculate_boundaries", bc_type)   
-       
-        select case(bc_type)
-        case("neumann")
-         
-            
-        case("dirichlet") 
-        !! Define the values at ice-ocean interface    
-        ! Get the surface_id of the ice-ocean interface
-            shape_option=option_shape(trim(melt_path)//"/melt_surfaceID")
-            allocate(surf_id(1:shape_option(1)))
-            call get_option(trim(melt_path)//'/melt_surfaceID',surf_id)
-            call allocate(surface_ids)
-            call insert(surface_ids,surf_id)
-            mesh => extract_mesh(state,"VelocityMesh")
-            ! Input, mesh, surface_id
-            ! Output surface_mesh,surface_element_list
-            call melt_surf_mesh(mesh,surface_ids,surface_mesh,surface_nodes,surface_element_list)
-
-            T => extract_scalar_field(state,"Temperature")
-            S => extract_scalar_field(state,"Salinity")
-            do i=1,size(surface_nodes)
-                the_node = surface_nodes(i)
-                
-                call set(T,the_node,0.0)
-                call set(S,the_node,34.0)
-               
-            enddo
-!            T_bc => extract_scalar_field(state,"Tb")
-!            S_bc => extract_scalar_field(state,"Sb")
-        case default
-            FLAbort('Unknown BC for TKE')
-        end select 
-    
-
-    ewrite(1,*) "---------End melt_surf_init---------------------------------"
-
- end subroutine melt_surf_init
-
-subroutine melt_allocate_surface(state)
-
-    type(state_type), intent(inout)     :: state
-    type(vector_field), pointer         :: positions   
-    type(mesh_type), pointer            :: mesh
-    type(mesh_type)                     :: surface_mesh
-    type(integer_set)                   :: surface_elements
-    integer, dimension(:), pointer      :: surface_nodes ! allocated and returned by create_surface_mesh
-    integer, dimension(:), allocatable  :: surface_element_list
-    integer    :: ele, face,i,j,k
-!! The local coordinates of the coordinate in the owning element
-!!    real, dimension(:), allocatable    :: local_coord, local_coord_surf
-    real, dimension(:), allocatable     :: coord
-! for transform_facet_to_physical
-    real, dimension(:,:), allocatable   :: normal
-    real, dimension(:), allocatable     :: av_normal !Average of normal
-!!    integer, dimension(:), pointer      :: ele_faces,face_neighs
-    type(element_type), pointer     :: x_shape_f
-    real                                :: al,be,ga
-    real, dimension(:), allocatable     :: xyz !New location of surface_mesh, dist_meltrate away from the boundary
-!integer, save         :: melt_surfaceID
-    character(len=OPTION_PATH_LEN)      :: path
-    type(integer_set)                   :: surface_ids
-    integer, dimension(:),allocatable   :: surf_id
-
-! For the vector averaging schem, taking the area of the surface element etc 
-    real, dimension(:,:), allocatable   :: table
-    integer, dimension(:,:), allocatable :: node_occupants
-    integer                             :: st,en,node,dim_vec
-    real                                :: area_sum
-    integer, dimension(:), allocatable  :: sf_nodes_ar
-    integer, dimension(2) :: shape_option
-
-    ewrite(1,*) "-------Begin melt_allocate_surface---------"
-    path = "/ocean_forcing/iceshelf_meltrate/Holland08"
-    ! Get the surface_id of the ice-ocean interface
-    shape_option=option_shape(trim(path)//"/melt_surfaceID")
-    allocate(surf_id(1:shape_option(1)))
-    call get_option(trim(path)//'/melt_surfaceID',surf_id)
-    call allocate(surface_ids)
-    call insert(surface_ids,surf_id)
-!    deallocate(surf_id)
-!    ewrite(1,*) "aft_surface_ids: ", set2vector(surface_ids)
-
-    mesh => extract_mesh(state,"VelocityMesh")
-! Input, mesh, surface_id
-! Output surface_mesh,surface_element_list
-    call melt_surf_mesh(mesh,surface_ids,surface_mesh,surface_nodes,surface_element_list)
-
-    positions => extract_vector_field(state,"Coordinate")
-!SINK THE surface_mesh!!
-!   call get_option("/geometry/dimension/", dimension)
-    call allocate(surface_positions,positions%dim,mesh,"MySurfacePosition")
-    call allocate(funky_positions,surface_positions%dim, surface_positions%mesh,"MyFunkyPosition")
-    
-
-! Remap the positions vector to the surface mesh, surface_positions
-    call remap_field_to_surface(positions, surface_positions, surface_element_list)  
-   
-! Loop over # of surface elements
-!nf = size(surface_element_list)
-!2 comes because there are 2 nodes per a surface_element (assumption!)
-! this number could be three for 3D problem, room for change
-! dim_vec is the dimension of the positions, either 2 or 3.
-! node_occupants(1,1:dim_vec*nf) = nodes that occupies the surface elements
-! node_occupants(2,:) = surface elements
-! table(1,:) = area of the face
-! table(2,:) = normal_x
-! table(3,:) = normal_y
-! table(4,:) = normal_z
-    dim_vec = positions%dim
-    allocate(coord(dim_vec))
-    allocate(node_occupants(2,dim_vec*size(surface_element_list)))
-    allocate(table(dim_vec+1,dim_vec*size(surface_element_list)))
-    allocate(av_normal(dim_vec))
-    allocate(xyz(dim_vec))
-   
-    do i=1,size(surface_element_list)        
-        st = 1+dim_vec*(i-1)
-        en = dim_vec+dim_vec*(i-1)
-        face = surface_element_list(i)
-        node_occupants(2,st:en) = face
-        node_occupants(1,st:en) = face_global_nodes(positions, face)
-        ! Calculate the area of the surface element
-        ! For 2D, the area is the length
-        if (dim_vec .eq. 2) then
-            table(1,st:en) = calc_area2(positions,node_occupants(1,st:en))
-        endif
-        
-        ! For 3D, the area become the area of the triangle (assumption here).
-        ! Subroutine calc_area3 uses Heron's formula.
-        if (dim_vec .eq. 3) then
-            table(1,st:en) = calc_area3(positions,node_occupants(1,st:en)) 
-        endif
-        ! Compute the normal vector at the surface element.
-        x_shape_f=>face_shape(positions,face)
-        allocate(normal(dim_vec,x_shape_f%ngi)) !(dim x x_shape_f%ngi)
-        call transform_facet_to_physical(positions, face, normal=normal)
-        av_normal = sum(normal,2) 
-        !"normal" should be the same; since the normal vector on the quadrature point does not change.
-        av_normal = av_normal/(sum(av_normal*av_normal))**(0.5) ! normalize the vector
-        deallocate(normal)
-        !Since av_normal is outward to the boundary, we will put -sign. We want the vector into the boundary
-        av_normal = -av_normal
-        ! Store av_normal, the normal vector averaged over quadrature points, in table        
-        do j=1,dim_vec
-            table(j+1,st:en)=av_normal(j)
-        enddo
-    enddo
-!! Now loop over surface_nodes   
-!        ewrite(1,*) "table(1,1:3): ", table(1,1:3)
-!        ewrite(1,*) "table(2,1:3),normal_x: ", table(2,1:3)
-!        ewrite(1,*) "table(3,1:3),normal_y: ", table(3,1:3)
-!        ewrite(1,*) "table(4,1:3),normal_z: ", table(4,1:3)
-!!In this loop, we will average adjacent normal vectors, using the area of the surface elements.
-    
-    !allocate(sf_nodes_ar(size(node_occupants(1,:))))
-    call allocate(sf_nodes)
-!    call insert(sf_nodes,sf_nodes_ar)
-    do i=1,size(node_occupants(1,:))
-        node = node_occupants(1,i)
-        av_normal(:) = 0.0
-        area_sum = 0.0
-        ! This loop average the normal vector based on the areas of surface elements
-        !, which share the "node". Using hash table may speed up this process
-        do j=1,size(node_occupants(1,:))
-            !Pick the surface elements that sheare the "node"
-            if (node_occupants(1,j) .eq. node) then 
-                do k=1,dim_vec
-                    !table(1,j) = area of the surface element
-                    av_normal(k) = av_normal(k) + table(k+1,j)*table(1,j)               
-                enddo 
-                !The total areas of the surface elements that occupies the "node"               
-                area_sum = area_sum + table(1,j)  
-            endif             
-        enddo
-    
-        av_normal = av_normal / area_sum
-        ! normalize 
-        av_normal = av_normal/(sum(av_normal*av_normal))**(0.5) 
-         
-        dist_meltrate = abs(dist_meltrate)
-        ! Shift the location of the surface nodes.
-        !The coordinate of the surface node.
-        
-        coord = node_val(positions,node) 
-        
-        ! dist_meltrate = ||xyz - coord|| xyz and coord are vectors
-        xyz = av_normal*dist_meltrate + coord
-        
-        call set(funky_positions,node,xyz) ! Set the coordinate of sinked nodes, funky positions.
-       
-        call set(surface_positions,node,coord) !Original coordinate of the surface
-!!        ewrite(1,*) "--------------------------------"
-!!        ewrite(1,*) "node: ", node
-!!        ewrite(1,*) "av_normal: ", av_normal
-!!        ewrite(1,*) "funky: ", xyz
-!!        ewrite(1,*) "coord: ", coord
-!!        ! Save the corresponding node number.
-!!        !sf_nodes_ar(i) = node
-        call insert(sf_nodes,node)
-       
-        node = 0
-    enddo
-     
-
-    deallocate(coord)
-    deallocate(table)
-    deallocate(node_occupants)    
-    deallocate(av_normal)
-    deallocate(xyz)
-    call deallocate(surface_ids)
-     ewrite(1,*) "-------End melt_allocate_surface---------"
-end subroutine melt_allocate_surface
-
-
-
-! Calculate the melt rate
- subroutine melt_surf_calc(state)
-
-    type(state_type), intent(inout)     :: state
-    type(scalar_field), pointer         :: Tb, Sb,MeltRate,Heat_flux,Salt_flux
-    type(scalar_field), pointer         :: scalarfield
-    type(vector_field), pointer         :: velocity,positions  
-    !Debugging purpose pointer
-    type(scalar_field), pointer         :: T_loc,S_loc,P_loc
-    type(vector_field), pointer         :: V_loc,Location,Location_org 
-    real, dimension(:), allocatable     :: vel
-    integer                             :: i,j,k
-    ! Some internal variables  
-    real                                :: speed, T,S,P,Aa,Bb,Cc,topo
-    real                                ::loc_Tb,loc_Sb,loc_meltrate,loc_heatflux,loc_saltflux
-    ! Aa*Sb^2+Bv*Sb+Cc
-    real :: arg = -1.0
-    !Sink mesh part
-    integer                             :: ele,face,node,stat,the_node
-    real, dimension(:), allocatable     :: local_coord,coord
-    type(element_type), pointer         :: x_shape_f
-    integer, dimension(:), allocatable  :: surface_node_list,node_lists 
-    type(scalar_field)                  :: re_temperature,re_salinity,re_pressure
-    type(vector_field)                  :: re_velocity
-  
-    
-    ewrite(1,*) "-------Begin melt_surf_calc------------"
-
-    MeltRate => extract_scalar_field(state,"MeltRate")
-    Tb => extract_scalar_field(state,"Tb")
-    Sb => extract_scalar_field(state,"Sb")
-    call set(MeltRate,setnan(arg))
-    call set(Tb,setnan(arg))
-    call set(Sb,setnan(arg))
-    Heat_flux => extract_scalar_field(state,"Heat_flux")
-    Salt_flux => extract_scalar_field(state,"Salt_flux")
-    call set(Heat_flux,setnan(arg))
-    call set(Salt_flux,setnan(arg))
-    
-! Debugging
-    T_loc => extract_scalar_field(state,"Tloc")
-    S_loc => extract_scalar_field(state,"Sloc")
-    P_loc => extract_scalar_field(state,"Ploc")
-    V_loc => extract_vector_field(state,"Vloc")
-    Location => extract_vector_field(state,"Location")
-    Location_org => extract_vector_field(state,"Location_org")
-    call set(T_loc,setnan(arg))
-    call set(S_loc,setnan(arg))
-    call set(P_loc,setnan(arg))
-    allocate(vel(V_loc%dim))
-    vel = setnan(arg)
-    call set(V_loc,vel)
-    call set(Location,vel)
-    call set(Location_org,vel)
-    
-    positions => extract_vector_field(state,"Coordinate")
-    !my positions
-  
-      ! Surface node list
-    allocate(surface_node_list(key_count(sf_nodes)))
-    ! Make it to vector from integer_set.
-    ! sf_nodes is calculated in "melt_allocate_surface"
-    surface_node_list=set2vector(sf_nodes)
-    
-    ! Remap temperature, salinity, pressure, and velocity onto positions mesh
-    scalarfield => extract_scalar_field(state,"Temperature")
-    
-    call allocate(re_temperature,positions%mesh,name="ReTemperature")
-    
-    call remap_field(scalarfield,re_temperature,stat)
-   
-    ! Salinity
-    scalarfield=> extract_scalar_field(state,"Salinity")
- 
-    call allocate(re_salinity,positions%mesh, name="ReSalinity")
-    call remap_field(scalarfield,re_salinity,stat)
-
-    ! Pressure
-    scalarfield => extract_scalar_field(state,"Pressure")
-    call allocate(re_pressure,positions%mesh, name="RePressure")
-    call remap_field(scalarfield,re_pressure,stat)
-   
-    ! Velocity    
-    velocity => extract_vector_field(state,"Velocity")
-    call allocate(re_velocity,velocity%dim,positions%mesh,name="ReVelocity")
-    call remap_field(velocity, re_velocity, stat)
- 
-    allocate(local_coord(positions%dim+1))
-    allocate(coord(positions%dim))
-   
-  
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!
-     !surface_positions
-
-    !! Loope over the surface nodes to calculate melt rate etc.
-    do i=1,size(surface_node_list)
-        the_node = surface_node_list(i)
-        !!! Interpolating   
-        coord = node_val(funky_positions,the_node)
-        call picker_inquire(positions, coord, ele, local_coord,global=.true.)
-
-        !! If sum(local_coord) is not equal to 1,      
-        !! we know that this coord (funky position) does not exist in the domain.
-        
-        if (sum(local_coord) .gt. 2.0) then 
-            ewrite(1,*) "Funk coord: ", node_val(funky_positions,the_node)
-            !! node_val(surface_positions,the_node) = node_val(positions,the_node)
-            ewrite(1,*) "Original ele: ",ele   
-            ewrite(1,*) "sum of local_coord: ",  sum(local_coord)    
-            FLExit("Your funky_positions is out of the domain. Change melt_LayerLength.")
-        endif
-        !Number of nodes per element       
-        allocate(node_lists(ele_and_faces_loc(positions, ele))) 
-        node_lists =  ele_nodes(positions,ele) !Lists of nodes that occupies the element.
-        ! This method of finding values at funky_positions works for P1, which are T,S, and velocity.
-        coord = 0.0
-        T = 0.0
-        S = 0.0
-        speed = 0.0
-        vel = 0.0
-
-        do j=1,size(node_lists)
-            node = node_lists(j)
-            coord = coord + node_val(positions,node)*local_coord(j)
-            T = T + node_val(re_temperature,node)*local_coord(j)
-            S = S + node_val(re_salinity,node)*local_coord(j)
-            vel = vel + node_val(re_velocity,node)*local_coord(j)
-        enddo
-        deallocate(node_lists)
-        ! Luckly P needs to be at the surface for the three equations
-        P = node_val(re_pressure,the_node)
-        speed = sqrt(sum(vel**2))
-        
-        if (speed .lt. 0.001) then
-            speed = 0.001
-             !ewrite(1,*) "----------iceshelf, speed less----", the_node,speed
-        endif
-        topo = -7.53e-8*P ! constant = -7.53e-8 [C Pa^(-1)] comes from Holland and Jenkins Table 1
-        
-        !! Define Aa,Bb,Cc
-        !! Aa*Sb**2 + Bb*Sb + Cc = 0.0
-        Aa = -gammaS*speed*cI*a + a*c0*gammaT*speed
-        
-        Bb = -gammaS*speed*L + gammaS*speed*S*cI*a
-        Bb = Bb - gammaS*speed*cI*(b+topo) + gammaS*speed*cI*TI
-        Bb = Bb - c0*gammaT*speed*T + c0*gammaT*speed*(b+topo)
-
-        Cc = gammaS*speed*S*L +gammaS*speed*S*cI*(b+topo) + gammaS*speed*S*(-cI*TI)
-        
-        
-        !! This could be a linear equation if Aa=0
-        if (Aa .eq. 0.0) then
-            loc_Sb = -Cc/Bb
-        else
-            !! Calculate for the 2nd oewrite(1,*) "size(surface_element_list)"rder polynomial.
-            !! We have two solutions. 
-            loc_Sb = (-Bb + sqrt(Bb**2 - 4.0*Aa*Cc))/(2.0*Aa)
-            !! loc_Sb has to be larger than 0; since the salinity in the ocean is positive definite.
-            if (loc_Sb .lt. 0.0) then
-                loc_Sb = (-Bb - sqrt(Bb**2 - 4.0*Aa*Cc))/(2.0*Aa)
-            endif
-        endif
-        !ewrite(1,*) "----------iceshelf, loc_Sb----", loc_Sb
-        loc_Tb = a*loc_Sb + b + topo
-        loc_meltrate = gammaS*speed*(S-loc_Sb)/loc_Sb
-        !! Heat flux to the ocean
-        loc_heatflux = c0*(gammaT*speed+ loc_meltrate)*(T-loc_Tb) ! or loc_meltrate*L + loc_meltrate*cI*(loc_Tb-TI)
-        loc_saltflux = (gammaS*speed+loc_meltrate)*(S-loc_Sb)
-        !! Some debugging
-        !!ewrite(1,*) "melt_rate: ",loc_meltrate
-        !!ewrite(1,*) "tLHS: ", c0*gammaT*speed*(T-loc_Tb)
-        !!ewrite(1,*) "tRHS: ", loc_meltrate*L + loc_meltrate*cI*(loc_Tb-TI)
-        !!ewrite(1,*) "sLHS: ", gammaS*speed*(S-loc_Sb)
-        !!ewrite(1,*) "sRHS: ", loc_meltrate*loc_Sb
-        !!ewrite(1,*) "bLHS: ", loc_Tb
-        !!ewrite(1,*) "bRHS: ", a*loc_Sb + b + topo
-
-        !! These are needed to implement BCs.
-        call set(MeltRate, the_node, loc_meltrate)
-        call set(Tb, the_node, loc_Tb)
-        call set(Sb, the_node, loc_Sb)
-        call set(Heat_flux, the_node,loc_heatflux)
-        call set(Salt_flux, the_node,loc_saltflux)
-        !! More or less for debugging purposes.
-        call set(T_loc,the_node,T)
-        call set(S_loc,the_node,S)
-        call set(P_loc,the_node,P)
-        call set(V_loc,the_node,vel)
-        call set(Location,the_node,node_val(funky_positions,the_node))
-        call set(Location_org,the_node,node_val(positions,the_node))
-        !! BC test
-        !!call set(TT,the_node,11.1)
-!         ewrite(1,*) "----------iceshelf, loc_saltflux----", the_node,loc_saltflux
-!        ewrite(1,*) "----------melt_surf_calc, loc_Tb----", the_node,loc_Tb
-!        ewrite(1,*) "----------iceshelf, surfaceS----",node_val(re_salinity,the_node)
-!        ewrite(1,*) "----------line 477 iceshelf, loc_meltrate----",loc_meltrate
-!        ewrite(1,*) "----------iceshelf, node_val(TT,the_node)----",node_val(TT,the_node)
-    enddo
-
-    deallocate(local_coord)
-    deallocate(coord)
-    ewrite(1,*) "-----END melt_surf_calc-------"
-end subroutine melt_surf_calc
-
-
-subroutine melt_bc(state)
-    type(state_type), intent(inout)     :: state
-  !! BC
-    type(scalar_field), pointer         :: TT,SS
-    type(scalar_field), pointer         :: scalar_surface
-    type(mesh_type), pointer            :: ice_mesh
-    character(len=FIELD_NAME_LEN)       :: bc_type
-    type(scalar_field), pointer         :: T_bc,S_bc
-    integer, dimension(:), allocatable  :: surface_node_list
-    integer                             :: i, the_node
-!! Insert BC for temperature and salinity. This could be a separate subroutine
-
-    
-        TT=> extract_scalar_field(state,"Temperature")
-        SS => extract_scalar_field(state,"Salinity")
-        
-        !This change with bc type
-        call get_option("/ocean_forcing/iceshelf_meltrate/Holland08/calculate_boundaries", bc_type)   
-       
-        select case(bc_type)
-        case("neumann")
-            T_bc => extract_scalar_field(state,"Heat_flux")
-            S_bc => extract_scalar_field(state,"Salt_flux")
-            do i=1,node_count(T_bc)
-                call set(T_bc,node_val(T_bc,i)*10.0**3)
-                call set(S_bc,node_val(S_bc,i)*10.0**3)
-            enddo
-            
-        case("dirichlet") 
-            T_bc => extract_scalar_field(state,"Tb")
-            S_bc => extract_scalar_field(state,"Sb")
-        case default
-            FLAbort('Unknown BC for TKE')
-        end select 
-    
-    
-        ! Surface node list
-        allocate(surface_node_list(key_count(sf_nodes)))
-        ! Make it to vector from integer_set.
-        ! sf_nodes is calculated in "melt_allocate_surface"
-        surface_node_list=set2vector(sf_nodes)
-        
-        ! create a surface mesh to place values onto. This is for the top surface
-        call get_boundary_condition(TT, 'temperature_iceshelf_BC', surface_mesh=ice_mesh) 
-        call allocate(ice_surfaceT, ice_mesh, name="ice_surfaceT")
-        call get_boundary_condition(SS, 'salinity_iceshelf_BC', surface_mesh=ice_mesh) 
-        call allocate(ice_surfaceS, ice_mesh, name="ice_surfaceS")
-        ! Define ice_surfaceT according to Heat_flux?
-        ewrite(1,*) "node_count(ice_surfaceT)",node_count(ice_surfaceT)
-        do i=1,node_count(ice_surfaceT)
-            the_node = surface_node_list(i)
-            call set(ice_surfaceT,i,node_val(T_bc,the_node))
-            call set(ice_surfaceS,i,node_val(S_bc,the_node))
-        enddo
-        !! Temperature
-        scalar_surface => extract_surface_field(TT, 'temperature_iceshelf_BC', "value")    
-        call remap_field(ice_surfaceT, scalar_surface)
-        !! Salinity
-        scalar_surface => extract_surface_field(SS, 'salinity_iceshelf_BC', "value")    
-        call remap_field(ice_surfaceS, scalar_surface)
-!        ewrite(1,*) "iceBC, ice_element_list",ice_element_list 
-!        ewrite(1,*) "iceBC, node_count(scalar_surface)", node_count(scalar_surface)
-!        ewrite(1,*) "node_val(scalar_surface,1): ", node_val(scalar_surface,1)
-
-!! Insert BC for salinity
-      
-
-end subroutine melt_bc
-
-
-
-!!------------------------------------------------------------------!
-!!                       Private subroutines                        !
-!!------------------------------------------------------------------!
-
-subroutine melt_surf_mesh(mesh,surface_ids,surface_mesh,surface_nodes,surface_element_list)
-    
-    type(mesh_type), intent(in)                       :: mesh
-    type(integer_set), intent(in)                     :: surface_ids
-    type(mesh_type), intent(out)                      :: surface_mesh
-    integer, dimension(:), pointer, intent(out)       :: surface_nodes ! allocated and returned by create_surface_mesh
-    integer, dimension(:), allocatable, intent(out)   :: surface_element_list
-    !    integer, dimension(:), allocatable                :: surface_nodes_out
-    type(scalar_field)                                :: surface_field   
-    type(integer_set)                                 :: surface_elements
-    integer :: i
-    ! create a set of surface elements that have surface id in 'surface_ids'
-
-    call allocate(surface_elements)
-    do i=1, surface_element_count(mesh)
-!        ewrite(1,*) "surf_normal surface_element_id(mesh, i)", surface_element_id(mesh, i)
-!        ewrite(1,*) "surf_normal surface_ids", set2vector(surface_ids)
-        if (has_value(surface_ids, surface_element_id(mesh, i))) then
-             call insert(surface_elements, i)
-             
-        end if
-    end do
-    
-    allocate(surface_element_list(key_count(surface_elements)))
-    surface_element_list=set2vector(surface_elements)
-    call create_surface_mesh(surface_mesh, surface_nodes, mesh, surface_element_list, name=trim(mesh%name)//"ToshisMesh")
-    
-
-end subroutine melt_surf_mesh
-
-
-real function setnan(arg)
-     real :: arg
-        setnan = sqrt(arg)
-end function
-
-real function calc_area2(positions,nodes)
-     type(vector_field), pointer        :: positions
-     integer, dimension(2)              :: nodes
-     real, dimension(2)                 :: coord1,coord2
-!     real                               :: area
-     coord1 = node_val(positions,nodes(1))
-     coord2 = node_val(positions,nodes(2))
-    
-     calc_area2 = (sum((coord2 - coord1)**2))**0.5
-     
-end function
-
-real function calc_area3(positions,nodes)
-     type(vector_field), pointer        :: positions
-     integer, dimension(3)              :: nodes
-     real, dimension(3)                 :: coord1,coord2,coord3
-     real                               :: a,b,c,s
-!     real                               :: area
-     coord1 = node_val(positions,nodes(1))
-     coord2 = node_val(positions,nodes(2))
-     coord3 = node_val(positions,nodes(3))
-! Use Heron's formula to calculate the area of triangle
-     a = (sum((coord2 - coord1)**2))**0.5
-     b = (sum((coord3 - coord1)**2))**0.5
-     c = (sum((coord2 - coord3)**2))**0.5
-     s = 0.5*(a+b+c)     
-     calc_area3 = (s*(s-a)*(s-b)*(s-c))**0.5
-
-end function
-
-end module iceshelf_meltrate_surf_normal

=== removed file 'schemas/multiphase_interaction.rnc'
--- schemas/multiphase_interaction.rnc	2012-06-15 19:32:03 +0000
+++ schemas/multiphase_interaction.rnc	1970-01-01 00:00:00 +0000
@@ -1,58 +0,0 @@
-# Phase interaction options for Fluidity's multiphase flow model.
-
-multiphase_interaction = 
-   (
-      ## Phase interaction options for Fluidity's multiphase flow model.
-      element multiphase_interaction {
-         comment,
-         fluid_particle_drag?,
-         heat_transfer?
-      }
-   ) 
-
-fluid_particle_drag =
-   (
-      ## Fluid-particle drag term.
-      element fluid_particle_drag {
-         drag_correlation
-      }
-   )
-
-drag_correlation =
-   (
-      ## Stokes drag correlation: 24/Re_p
-      ## where Re_p is the particle Reynolds number.
-      element drag_correlation {
-         attribute name { "stokes" },
-         comment
-      }|
-      ## Fluid-particle drag term by Wen & Yu (1966).
-      element drag_correlation {
-         attribute name { "wen_yu" },
-         comment
-      }|
-      ## Fluid-particle drag term by Ergun (1952).
-      element drag_correlation {
-         attribute name { "ergun" },
-         comment
-      }
-   )
-   
-heat_transfer =
-   (
-      ## Heat transfer term for the 
-      ## multiphase internal energy equation.
-      ## Note: Only for fluid-particle phase pairs.
-      element heat_transfer {
-         heat_transfer_coefficient
-      }
-   )
-
-heat_transfer_coefficient =
-   (
-      ## Heat transfer coefficient by Gunn (1978).
-      element heat_transfer_coefficient {
-         attribute name { "gunn" },
-         comment
-      }
-   )
\ No newline at end of file

=== removed file 'schemas/multiphase_interaction.rng'
--- schemas/multiphase_interaction.rng	2012-06-15 19:32:03 +0000
+++ schemas/multiphase_interaction.rng	1970-01-01 00:00:00 +0000
@@ -1,65 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<grammar xmlns:a="http://relaxng.org/ns/compatibility/annotations/1.0" xmlns="http://relaxng.org/ns/structure/1.0">
-  <!-- Phase interaction options for Fluidity's multiphase flow model. -->
-  <define name="multiphase_interaction">
-    <element name="multiphase_interaction">
-      <a:documentation>Phase interaction options for Fluidity's multiphase flow model.</a:documentation>
-      <ref name="comment"/>
-      <optional>
-        <ref name="fluid_particle_drag"/>
-      </optional>
-      <optional>
-        <ref name="heat_transfer"/>
-      </optional>
-    </element>
-  </define>
-  <define name="fluid_particle_drag">
-    <element name="fluid_particle_drag">
-      <a:documentation>Fluid-particle drag term.</a:documentation>
-      <ref name="drag_correlation"/>
-    </element>
-  </define>
-  <define name="drag_correlation">
-    <choice>
-      <element name="drag_correlation">
-        <a:documentation>Stokes drag correlation: 24/Re_p
-where Re_p is the particle Reynolds number.</a:documentation>
-        <attribute name="name">
-          <value>stokes</value>
-        </attribute>
-        <ref name="comment"/>
-      </element>
-      <element name="drag_correlation">
-        <a:documentation>Fluid-particle drag term by Wen &amp; Yu (1966).</a:documentation>
-        <attribute name="name">
-          <value>wen_yu</value>
-        </attribute>
-        <ref name="comment"/>
-      </element>
-      <element name="drag_correlation">
-        <a:documentation>Fluid-particle drag term by Ergun (1952).</a:documentation>
-        <attribute name="name">
-          <value>ergun</value>
-        </attribute>
-        <ref name="comment"/>
-      </element>
-    </choice>
-  </define>
-  <define name="heat_transfer">
-    <element name="heat_transfer">
-      <a:documentation>Heat transfer term for the 
-multiphase internal energy equation.
-Note: Only for fluid-particle phase pairs.</a:documentation>
-      <ref name="heat_transfer_coefficient"/>
-    </element>
-  </define>
-  <define name="heat_transfer_coefficient">
-    <element name="heat_transfer_coefficient">
-      <a:documentation>Heat transfer coefficient by Gunn (1978).</a:documentation>
-      <attribute name="name">
-        <value>gunn</value>
-      </attribute>
-      <ref name="comment"/>
-    </element>
-  </define>
-</grammar>

=== modified file 'schemas/prognostic_field_options.rnc'
--- schemas/prognostic_field_options.rnc	2012-08-22 16:03:09 +0000
+++ schemas/prognostic_field_options.rnc	2012-09-04 15:00:37 +0000
@@ -91,8 +91,7 @@
       ##  T the scalar field value on the surface,
       ##  C0 is the input order zero coefficient and
       ##  C1 is the input order one coefficient.
-      ##  THIS WILL ONLY WORK FOR CONTINUOUS GALERKIN OR CONTROL VOLUME 
-      ##  (USING AN ELEMENTGRADIENT DIFFUSION SCHEME) SPATIAL DISCRETISATIONS.
+      ##  THIS WILL ONLY WORK FOR CONTINUOUS GALERKIN SPATIAL DISCRETISATION
       element type {
          attribute name { "robin" },
          ##  The order zero coefficient represented as C0 in

=== modified file 'schemas/prognostic_field_options.rng'
--- schemas/prognostic_field_options.rng	2012-08-22 16:03:09 +0000
+++ schemas/prognostic_field_options.rng	2012-09-04 15:00:37 +0000
@@ -113,8 +113,7 @@
  T the scalar field value on the surface,
  C0 is the input order zero coefficient and
  C1 is the input order one coefficient.
- THIS WILL ONLY WORK FOR CONTINUOUS GALERKIN OR CONTROL VOLUME 
- (USING AN ELEMENTGRADIENT DIFFUSION SCHEME) SPATIAL DISCRETISATIONS.</a:documentation>
+ THIS WILL ONLY WORK FOR CONTINUOUS GALERKIN SPATIAL DISCRETISATION</a:documentation>
       <attribute name="name">
         <value>robin</value>
       </attribute>

=== modified file 'schemas/spatial_discretisation.rnc'
--- schemas/spatial_discretisation.rnc	2012-06-02 12:15:06 +0000
+++ schemas/spatial_discretisation.rnc	2012-09-04 15:00:37 +0000
@@ -494,7 +494,7 @@
             ## basis functions of the parent finite element mesh to
             ## form the divergence.
             ##
-            ## DOES NOT CURRENTLY WORK WITH WEAK DIRICHLET BOUNDARY CONDITIONS!
+            ## DOES NOT CURRENTLY WORK WITH ROBIN OR WEAK DIRICHLET BOUNDARY CONDITIONS!
             ##
             ## Based on schemes in Handbook of Numerical Analysis,
             ## P.G. Ciarlet, J.L. Lions eds, vol 7, pp 713-1020

=== modified file 'schemas/spatial_discretisation.rng'
--- schemas/spatial_discretisation.rng	2012-06-02 12:15:06 +0000
+++ schemas/spatial_discretisation.rng	2012-09-04 15:00:37 +0000
@@ -546,7 +546,7 @@
 basis functions of the parent finite element mesh to
 form the divergence.
 
-DOES NOT CURRENTLY WORK WITH WEAK DIRICHLET BOUNDARY CONDITIONS!
+DOES NOT CURRENTLY WORK WITH ROBIN OR WEAK DIRICHLET BOUNDARY CONDITIONS!
 
 Based on schemes in Handbook of Numerical Analysis,
 P.G. Ciarlet, J.L. Lions eds, vol 7, pp 713-1020</a:documentation>

=== modified file 'tests/diff_src_robin_1d/diff_src_robin_1d.xml'
--- tests/diff_src_robin_1d/diff_src_robin_1d.xml	2012-06-02 12:15:06 +0000
+++ tests/diff_src_robin_1d/diff_src_robin_1d.xml	2012-09-04 15:00:37 +0000
@@ -6,8 +6,7 @@
   <tags>flml</tags>
   <problem_definition length="short" nprocs="1">
     <command_line>
-../../bin/fluidity diff_src_robin_1d_cg_p1_A.flml
-../../bin/fluidity diff_src_robin_1d_cvelegrad_p1_A.flml 
+../../bin/fluidity diff_src_robin_1d_cg_p1_A.flml 
     </command_line>
     <!-- One dimensional problem for diffusion, source with right robin bc and left neumann zero compared to analytic-->
   </problem_definition>
@@ -20,14 +19,6 @@
 from fluidity_tools import stat_parser as stat
 cg_p1_A_error_linf = stat("diff_src_robin_1d_cg_p1_A.stat")["fluid"]["TError"]["max"][-1]
     </variable>
-    <variable name="cvelegrad_p1_A_error_l2" language="python">
-from fluidity_tools import stat_parser as stat
-cvelegrad_p1_A_error_l2 = stat("diff_src_robin_1d_cvelegrad_p1_A.stat")["fluid"]["TError"]["cv_l2norm"][-1]
-    </variable>
-    <variable name="cvelegrad_p1_A_error_linf" language="python">
-from fluidity_tools import stat_parser as stat
-cvelegrad_p1_A_error_linf = stat("diff_src_robin_1d_cvelegrad_p1_A.stat")["fluid"]["TError"]["max"][-1]
-    </variable>
   </variables>
   <pass_tests>
     <test name="Check cg_p1_A_error_l2 less than 1.0e-8" language="python">
@@ -36,12 +27,6 @@
     <test name="Check cg_p1_A_error_linf less than 1.0e-8" language="python">
 assert(abs(cg_p1_A_error_linf) &lt; 1.0e-8)
     </test>
-    <test name="Check cvelegrad_p1_A_error_l2 less than 1.0e-8" language="python">
-assert(abs(cvelegrad_p1_A_error_l2) &lt; 1.0e-8)
-    </test>
-    <test name="Check cvelegrad_p1_A_error_linf less than 1.0e-8" language="python">
-assert(abs(cvelegrad_p1_A_error_linf) &lt; 1.0e-8)
-    </test>
   </pass_tests>
   <warn_tests>
   </warn_tests>

=== removed file 'tests/diff_src_robin_1d/diff_src_robin_1d_cvelegrad_p1_A.flml'
--- tests/diff_src_robin_1d/diff_src_robin_1d_cvelegrad_p1_A.flml	2012-06-02 12:18:06 +0000
+++ tests/diff_src_robin_1d/diff_src_robin_1d_cvelegrad_p1_A.flml	1970-01-01 00:00:00 +0000
@@ -1,244 +0,0 @@
-<?xml version='1.0' encoding='utf-8'?>
-<fluidity_options>
-  <simulation_name>
-    <string_value lines="1">diff_src_robin_1d_cvelegrad_p1_A</string_value>
-  </simulation_name>
-  <problem_type>
-    <string_value lines="1">fluids</string_value>
-  </problem_type>
-  <geometry>
-    <dimension>
-      <integer_value rank="0">1</integer_value>
-    </dimension>
-    <mesh name="CoordinateMesh">
-      <from_file file_name="diff_src_robin_1d_A">
-        <format name="triangle"/>
-        <stat>
-          <include_in_stat/>
-        </stat>
-      </from_file>
-    </mesh>
-    <quadrature>
-      <degree>
-        <integer_value rank="0">2</integer_value>
-      </degree>
-    </quadrature>
-  </geometry>
-  <io>
-    <dump_format>
-      <string_value>vtk</string_value>
-    </dump_format>
-    <dump_period_in_timesteps>
-      <constant>
-        <integer_value rank="0">0</integer_value>
-      </constant>
-    </dump_period_in_timesteps>
-    <output_mesh name="CoordinateMesh"/>
-    <stat/>
-  </io>
-  <timestepping>
-    <current_time>
-      <real_value rank="0">0.0</real_value>
-    </current_time>
-    <timestep>
-      <real_value rank="0">2.5</real_value>
-      <comment>Time stepping is not actually needed for this test case but this is just to check the acceleration field solver</comment>
-    </timestep>
-    <finish_time>
-      <real_value rank="0">10.0</real_value>
-      <comment>Time stepping is not actually needed for this test case but this is just to check the acceleration field solver</comment>
-    </finish_time>
-  </timestepping>
-  <physical_parameters/>
-  <material_phase name="fluid">
-    <vector_field name="Velocity" rank="1">
-      <prescribed>
-        <mesh name="CoordinateMesh"/>
-        <value name="WholeMesh">
-          <constant>
-            <real_value shape="1" dim1="dim" rank="1">0.0</real_value>
-          </constant>
-        </value>
-        <output/>
-        <stat>
-          <include_in_stat/>
-        </stat>
-        <detectors>
-          <exclude_from_detectors/>
-        </detectors>
-      </prescribed>
-    </vector_field>
-    <scalar_field name="T" rank="0">
-      <prognostic>
-        <mesh name="CoordinateMesh"/>
-        <equation name="AdvectionDiffusion"/>
-        <spatial_discretisation>
-          <control_volumes>
-            <mass_terms>
-              <exclude_mass_terms/>
-            </mass_terms>
-            <face_value name="FirstOrderUpwind"/>
-            <diffusion_scheme name="ElementGradient"/>
-          </control_volumes>
-          <conservative_advection>
-            <real_value rank="0">1.0</real_value>
-          </conservative_advection>
-        </spatial_discretisation>
-        <temporal_discretisation>
-          <theta>
-            <real_value rank="0">1.0</real_value>
-          </theta>
-          <control_volumes>
-            <pivot_theta>
-              <real_value rank="0">1.0</real_value>
-            </pivot_theta>
-          </control_volumes>
-        </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-10</real_value>
-          </relative_error>
-          <absolute_error>
-            <real_value rank="0">1.0e-10</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">9.87654321</real_value>
-          </constant>
-        </initial_condition>
-        <boundary_conditions name="right">
-          <surface_ids>
-            <integer_value shape="1" rank="1">2</integer_value>
-          </surface_ids>
-          <type name="robin">
-            <order_zero_coefficient>
-              <constant>
-                <real_value rank="0">3.14</real_value>
-              </constant>
-            </order_zero_coefficient>
-            <order_one_coefficient>
-              <constant>
-                <real_value rank="0">0.78</real_value>
-              </constant>
-            </order_one_coefficient>
-          </type>
-        </boundary_conditions>
-        <tensor_field name="Diffusivity" rank="2">
-          <prescribed>
-            <value name="WholeMesh">
-              <isotropic>
-                <constant>
-                  <real_value rank="0">1.0</real_value>
-                </constant>
-              </isotropic>
-            </value>
-            <output/>
-          </prescribed>
-        </tensor_field>
-        <scalar_field name="Source" rank="0">
-          <prescribed>
-            <value name="WholeMesh">
-              <constant>
-                <real_value rank="0">1.234</real_value>
-              </constant>
-            </value>
-            <output/>
-            <stat/>
-            <detectors>
-              <exclude_from_detectors/>
-            </detectors>
-          </prescribed>
-        </scalar_field>
-        <output/>
-        <stat>
-          <include_cv_stats/>
-        </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="TAnalytic" rank="0">
-      <prescribed>
-        <mesh name="CoordinateMesh"/>
-        <value name="WholeMesh">
-          <python>
-            <string_value lines="20" type="code" language="python">def val(X,t):
-   # get the 1d coord
-   x = X[0]
-   
-   # define the analytic input
-   x_left      = 0.0
-   x_right     = 10.0
-   src         = 1.234
-   bc_right_r0 = 3.14
-   bc_right_r1 = 0.78  
-   
-   T_inf = bc_right_r0 / bc_right_r1   
-   gamma = bc_right_r1
-   
-   # calc the analytic solution
-   T = (src/2.0)*(x_right**2 - x**2) \
-     + (src*x_left)*(x - x_right) \
-     + (src/gamma)*(x_right - x_left) \
-     + T_inf 
-   
-   return T</string_value>
-          </python>
-        </value>
-        <output/>
-        <stat>
-          <include_cv_stats/>
-        </stat>
-        <detectors>
-          <exclude_from_detectors/>
-        </detectors>
-      </prescribed>
-    </scalar_field>
-    <scalar_field name="TError" rank="0">
-      <diagnostic>
-        <algorithm name="scalar_python_diagnostic" material_phase_support="single">
-          <string_value lines="20" type="code" language="python">field1=state.scalar_fields["T"]
-field2=state.scalar_fields["TAnalytic"]
-for i in range(field.node_count):
-   field.set(i, abs(field1.node_val(i)-field2.node_val(i)))</string_value>
-        </algorithm>
-        <mesh name="CoordinateMesh"/>
-        <output/>
-        <stat>
-          <include_cv_stats/>
-        </stat>
-        <convergence>
-          <include_in_convergence/>
-        </convergence>
-        <detectors>
-          <include_in_detectors/>
-        </detectors>
-        <steady_state>
-          <include_in_steady_state/>
-        </steady_state>
-      </diagnostic>
-    </scalar_field>
-  </material_phase>
-</fluidity_options>

=== removed directory 'tests/flux_bc'
=== removed file 'tests/flux_bc/Makefile'
--- tests/flux_bc/Makefile	2012-07-09 23:48:46 +0000
+++ tests/flux_bc/Makefile	1970-01-01 00:00:00 +0000
@@ -1,22 +0,0 @@
-preprocess:
-	@echo **********Creating meshes
-	gmsh -2 -o flux_bc_2d.msh src/flux_bc_2d.geo
-	gmsh -3 -o flux_bc_3d.msh src/flux_bc_3d.geo
-	../../bin/gmsh2triangle --2d flux_bc_2d.msh
-	../../bin/gmsh2triangle flux_bc_3d.msh
-
-run:
-	@echo **********Running 2D simulation
-	../../bin/fluidity flux_bc_2d.flml 
-	@echo **********Running 3D simulation
-	../../bin/fluidity flux_bc_3d.flml
-
-input: clean preprocess
-
-clean:
-	rm -f *.stat *.steady_state*
-	rm -f *.d.* *.vtu
-	rm -f *.msh
-	rm -f *.ele *.edge *.node *.poly *.face
-	rm -f matrixdump* *checkpoint* *.log* *.err*
-

=== removed file 'tests/flux_bc/flux_bc.xml'
--- tests/flux_bc/flux_bc.xml	2012-08-17 19:58:54 +0000
+++ tests/flux_bc/flux_bc.xml	1970-01-01 00:00:00 +0000
@@ -1,54 +0,0 @@
-<?xml version="1.0" encoding="UTF-8" ?>
-<!DOCTYPE testproblem SYSTEM "regressiontest.dtd">
-
-<testproblem>
-   
-  <name>flux_bc</name>
-  <owner userid="ctj10"/>
-  <tags>flml</tags>
-  
-  <problem_definition length="medium" nprocs="1">
-    <command_line>make run</command_line>
-    <!-- This checks the integral of the PhaseVolumeFraction field at t = 1 second. -->
-    <!-- The volumetric flux (r) is set to 0.5 m^3/s/m^2, and the flux boundary condition -->
-    <!-- enforces d/dt \int_{volume} vfrac = \int_{surface} r -->
-  </problem_definition>
-  
-  <variables>
-    <variable name="vfrac_integral_2d" language="python">
-from fluidity_tools import stat_parser
-s = stat_parser("flux_bc_2d.stat")
-vfrac_integral_2d = s["Tephra"]["PhaseVolumeFraction"]["integral"]
-    </variable>
-    
-    <variable name="vfrac_integral_3d" language="python">
-from fluidity_tools import stat_parser