Merge lp:~fluidity-core/fluidity/shallow-water-dev into lp:fluidity

[Dr Colin J Cotter]

This shallow-water code includes the changes necessary to correctly interpolate vector fields from the sphere for BDFM1 elements, plus a suite of test cases.

[Cian Wilson (cwilson)]

Hi Colin,

Only had a quick glance but can I check if these changes provide any new functionality to fluidity itself? I ask because you've changed mesh_options.rn[cg] and these options are now appearing in the fluidity schema too but I can't see how they're applicable to fluidity as I can only trace them back to the hybridized helmholtz solver. If they're not used in fluidity then I don't think they should appear in a piece of the schema that is visible within an flml but I may be missing something.

Unrelated to this merge but the parent option (geometry/mesh/from_mesh/constrain_type) for these new options is undocumented in the schema.

Cheers,
Cian

review: Needs Information

[Dr Colin J Cotter (colin-cotter)]

Hi Cian,
     These options are read in populate_state and make things happen in
Elements.F90 and various other places.

--cjc

On 09/11/11 16:17, Cian Wilson wrote:
> Review: Needs Information
>
> Hi Colin,
>
> Only had a quick glance but can I check if these changes provide any new functionality to fluidity itself? I ask because you've changed mesh_options.rn[cg] and these options are now appearing in the fluidity schema too but I can't see how they're applicable to fluidity as I can only trace them back to the hybridized helmholtz solver. If they're not used in fluidity then I don't think they should appear in a piece of the schema that is visible within an flml but I may be missing something.
>
> Unrelated to this merge but the parent option (geometry/mesh/from_mesh/constrain_type) for these new options is undocumented in the schema.
>
> Cheers,
> Cian

[lp:~fluidity-core/fluidity/shallow-water-dev]

3579. By Dr Colin J Cotter on 2011-11-09

A few fixes for adjoint, plus tweaking some swmls.

[Cian Wilson (cwilson)]

Hi Colin,

Could you give me line numbers in fluidity code where this happens please?

When I open up diamond to set up a new flml I get option of setting:
/geometry/mesh/from_mesh/constraint_type/check_continuity
and
/geometry/mesh/from_mesh/constraint_type/check_continuity_matrix.

However, when I grep for where these options are used I get:

cwilson@uisce:shallow-water-dev$ grep -ir check_continuity * -l
assemble/Hybridized_Helmholtz.F90
schemas/mesh_options.rnc
schemas/mesh_options.rng
tests/bdfm1_sphere/linear_williamson2_projection.swml

None of which appear to be fluidity related to me so shouldn't be available in an flml.

Cheers,
Cian

> Hi Cian,
> These options are read in populate_state and make things happen in
> Elements.F90 and various other places.
>
> --cjc
>
> On 09/11/11 16:17, Cian Wilson wrote:
> > Review: Needs Information
> >
> > Hi Colin,
> >
> > Only had a quick glance but can I check if these changes provide any new
> functionality to fluidity itself? I ask because you've changed
> mesh_options.rn[cg] and these options are now appearing in the fluidity schema
> too but I can't see how they're applicable to fluidity as I can only trace
> them back to the hybridized helmholtz solver. If they're not used in fluidity
> then I don't think they should appear in a piece of the schema that is visible
> within an flml but I may be missing something.
> >
> > Unrelated to this merge but the parent option
> (geometry/mesh/from_mesh/constrain_type) for these new options is undocumented
> in the schema.
> >
> > Cheers,
> > Cian

review: Needs Information

[Cian Wilson (cwilson)]

I can see where the already existing parent option /geometry/mesh/from_mesh/constraint_type is used in Populate_State just not the child options you are adding with this merge.

The parent option does need some documentation in the schema however.

[lp:~fluidity-core/fluidity/shallow-water-dev]

3580. By Patrick Farrell on 2011-11-10

Some XML fixes for Colin.

3581. By Patrick Farrell on 2011-11-10

Let's not put in a NonlinearVelocity in the adjoint run, for now.

3582. By Dr Colin J Cotter on 2011-11-11

Change place where option theta is looked for.

3583. By Dr Colin J Cotter on 2011-11-11

theta set.

3584. By Dr Colin J Cotter on 2011-11-11

Removing some dead tests and fixing some swmls.

[Dr Colin J Cotter (colin-cotter)]

Oh I see, fair point. I'll fix that.

all the best
--cjc

On 09/11/11 19:01, Cian Wilson wrote:
> I can see where the already existing parent option /geometry/mesh/from_mesh/constraint_type is used in Populate_State just not the child options you are adding with this merge.
>
> The parent option does need some documentation in the schema however.

[lp:~fluidity-core/fluidity/shallow-water-dev]

3585. By Dr Colin J Cotter on 2011-11-14

Move continuity debugging options to shallow_water_options.rnc as
requested by Cian.

3586. By Dr Colin J Cotter on 2011-11-14

Merge from trunk.

3587. By Dr Colin J Cotter on 2012-01-11

Merged from the trunk.

3588. By Dr Colin J Cotter on 2012-01-11

Forgot that this file has been renamed.

[David Ham (david-ham)]

I think this looks OK. Make sure it is current with trunk and check that all the major functionality is tested.

[lp:~fluidity-core/fluidity/shallow-water-dev]

3589. By Dr Colin J Cotter on 2012-02-13

Merged from trunk.

3590. By Dr Colin J Cotter on 2012-02-14

Uncommented code that I previously forgot to uncomment.

3591. By Dr Colin J Cotter on 2012-02-15

Get theta from correct place in adjoint callbacks.

3592. By Dr Colin J Cotter on 2012-02-15

Get theta from correct place in adjoint control callbacks too.

3593. By Dr Colin J Cotter on 2012-02-15

Some dead tests.

3594. By Dr Colin J Cotter on 2012-02-16

Remove pointer status from scalar field.

3595. By Dr Colin J Cotter on 2012-02-16

Don't output to the statfile twice at the end of the simulation.

This fixes shallow_water_adjoint_func_eval.

3596. By Dr Colin J Cotter on 2012-02-16

Merge from trunk.

3597. By Dr Colin J Cotter on 2012-02-21

Increased timesteps to 2 in test. This is not the solution, but I
can't work out why it doesn't work otherwise.

3598. By Dr Colin J Cotter on 2012-02-21

Need to calculate functional values before writing diagnostics.

3599. By Dr Colin J Cotter on 2012-02-27

Made makefiles and got a merge from the trunk.

3600. By Dr Colin J Cotter on 2012-02-28

Merge from trunk.

3601. By Dr Colin J Cotter on 2012-02-28

Error in merging makefiles.

3602. By Jemma Shipton on 2012-10-23

merging changes from Colin's bdfm1-nonlinear-sw branch

3603. By Jemma Shipton on 2012-10-23

get global numbering fix for quads from trunk

Unmerged revisions

3603. By Jemma Shipton on 2012-10-23

get global numbering fix for quads from trunk

3602. By Jemma Shipton on 2012-10-23

merging changes from Colin's bdfm1-nonlinear-sw branch

3601. By Dr Colin J Cotter on 2012-02-28

Error in merging makefiles.

3600. By Dr Colin J Cotter on 2012-02-28

Merge from trunk.

3599. By Dr Colin J Cotter on 2012-02-27

Made makefiles and got a merge from the trunk.

3598. By Dr Colin J Cotter on 2012-02-21

Need to calculate functional values before writing diagnostics.

3597. By Dr Colin J Cotter on 2012-02-21

Increased timesteps to 2 in test. This is not the solution, but I
can't work out why it doesn't work otherwise.

3596. By Dr Colin J Cotter on 2012-02-16

Merge from trunk.

3595. By Dr Colin J Cotter on 2012-02-16

Don't output to the statfile twice at the end of the simulation.

This fixes shallow_water_adjoint_func_eval.

3594. By Dr Colin J Cotter on 2012-02-16

Remove pointer status from scalar field.

=== modified file 'Makefile.in'
--- Makefile.in	2012-09-19 11:09:17 +0000
+++ Makefile.in	2012-10-23 15:04:25 +0000
@@ -186,7 +186,9 @@
 	@echo "    LD $(PROGRAM)"
 	@$(EVAL) $(LINKER) -o $(PROGRAM) main.o $(LIBS)
 
-sw: bin/shallow_water
+sw: bin/shallow_water 
+swm: bin/shallow_water_mimetic
+test_sw: bin/test_sw_advection
 be: bin/burgers_equation
 hh: bin/hybridized_helmholtz_solver
 
@@ -199,6 +201,22 @@
 	@$(EVAL) $(FLLINKER) -o bin/hybridized_helmholtz_solver main/Hybridized_Helmholtz_Solver.o $(LIBS)
 	@$(EVAL) $(FLLINKER) -o bin/shallow_water main/Shallow_Water.o $(LIBS)
 
+bin/shallow_water_mimetic: fluidity_library main/Shallow_Water_Mimetic.F90
+	@cd main; $(MAKE) Shallow_Water_Mimetic.o
+	@echo "BUILD shallow_water_mimetic"
+	@echo "    MKDIR bin"
+	@mkdir -p bin
+	@echo "    LD shallow_water_mimetic"
+	@$(EVAL) $(FLLINKER) -o bin/shallow_water_mimetic main/Shallow_Water_Mimetic.o $(LIBS)
+
+bin/test_sw_advection: fluidity_library main/test_sw_advection.F90
+	@cd main; $(MAKE) test_sw_advection.o
+	@echo "BUILD test_sw_advection"
+	@echo "    MKDIR bin"
+	@mkdir -p bin
+	@echo "    LD test_sw_advection"
+	@$(EVAL) $(LINKER) -o bin/test_sw_advection main/test_sw_advection.o $(LIBS)
+
 bin/hybridized_helmholtz_solver: fluidity_library main/Hybridized_Helmholtz_Solver.F90
 	@cd main; $(MAKE) Hybridized_Helmholtz_Solver.o
 	@echo "BUILD hybridized_helmholtz_solver"

=== modified file 'adjoint/Burgers_Control_Callbacks.F90'
--- adjoint/Burgers_Control_Callbacks.F90	2011-07-07 16:30:31 +0000
+++ adjoint/Burgers_Control_Callbacks.F90	2012-10-23 15:04:25 +0000
@@ -34,7 +34,7 @@
   use state_module
   use python_state
   use sparse_matrices_fields
-  use manifold_projections
+  use manifold_tools
   use mangle_dirichlet_rows_module
 
     implicit none

=== modified file 'adjoint/Makefile.dependencies'
--- adjoint/Makefile.dependencies	2011-07-19 16:55:03 +0000
+++ adjoint/Makefile.dependencies	2012-10-23 15:04:25 +0000
@@ -70,10 +70,9 @@
 Burgers_Control_Callbacks.o ../include/burgers_adjoint_controls.mod: \
    Burgers_Control_Callbacks.F90 ../include/fdebug.h ../include/fields.mod \
    ../include/global_parameters.mod \
-   ../include/mangle_dirichlet_rows_module.mod \
-   ../include/manifold_projections.mod ../include/python_state.mod \
-   ../include/sparse_matrices_fields.mod ../include/spud.mod \
-   ../include/state_module.mod
+   ../include/mangle_dirichlet_rows_module.mod ../include/manifold_tools.mod \
+   ../include/python_state.mod ../include/sparse_matrices_fields.mod \
+   ../include/spud.mod ../include/state_module.mod
 
 ../include/forward_main_loop.mod: Forward_Main_Loop.o
 	@true
@@ -128,7 +127,7 @@
    Shallow_Water_Adjoint_Callbacks.F90 ../include/adjoint_global_variables.mod \
    ../include/fdebug.h ../include/fields.mod ../include/global_parameters.mod \
    ../include/libadjoint_data_callbacks.mod ../include/mangle_options_tree.mod \
-   ../include/manifold_projections.mod ../include/populate_state_module.mod \
+   ../include/manifold_tools.mod ../include/populate_state_module.mod \
    ../include/sparse_matrices_fields.mod ../include/spud.mod \
    ../include/state_module.mod
 
@@ -140,7 +139,7 @@
    ../include/shallow_water_adjoint_controls.mod: \
    Shallow_Water_Control_Callbacks.F90 ../include/fdebug.h \
    ../include/fields.mod ../include/global_parameters.mod \
-   ../include/manifold_projections.mod ../include/python_state.mod \
+   ../include/manifold_tools.mod ../include/python_state.mod \
    ../include/sparse_matrices_fields.mod ../include/spud.mod \
    ../include/state_module.mod
 

=== modified file 'adjoint/Shallow_Water_Adjoint_Callbacks.F90'
--- adjoint/Shallow_Water_Adjoint_Callbacks.F90	2011-06-09 09:43:33 +0000
+++ adjoint/Shallow_Water_Adjoint_Callbacks.F90	2012-10-23 15:04:25 +0000
@@ -39,7 +39,7 @@
     use spud, only: get_option
     use adjoint_global_variables, only: adj_path_lookup
     use mangle_options_tree, only: adjoint_field_path
-    use manifold_projections
+    use manifold_tools
     use global_parameters, only: running_adjoint
     use populate_state_module
     implicit none
@@ -488,7 +488,7 @@
         path = adjoint_field_path(path)
       end if
 
-      call get_option(trim(path) // "/prognostic/temporal_discretisation/theta", theta)
+      call get_option("/timestepping/theta", theta)
       call get_option(trim(path) // "/prognostic/mean_layer_thickness", d0)
 
       if (hermitian==ADJ_FALSE) then
@@ -858,7 +858,7 @@
 
         ierr = adj_dict_find(adj_path_lookup, "Fluid::LayerThickness", path)
         call adj_chkierr(ierr)
-        call get_option(trim(path) // "/prognostic/temporal_discretisation/theta", theta)
+        call get_option("/timestepping/theta",theta)
         call get_option(trim(path) // "/prognostic/mean_layer_thickness", d0)
 
         eta_mesh => extract_mesh(matrices, "LayerThicknessMesh")

=== modified file 'adjoint/Shallow_Water_Control_Callbacks.F90'
--- adjoint/Shallow_Water_Control_Callbacks.F90	2011-06-11 13:26:55 +0000
+++ adjoint/Shallow_Water_Control_Callbacks.F90	2012-10-23 15:04:25 +0000
@@ -34,7 +34,7 @@
   use state_module
   use python_state
   use sparse_matrices_fields
-  use manifold_projections
+  use manifold_tools
 
     implicit none
 
@@ -154,7 +154,7 @@
                 adj_vfield => extract_vector_field(states(state_id), "AdjointLocalVelocityDelta") 
                 ! We need the AdjointVelocity for the option path only
                 adj_sfield2 => extract_scalar_field(states(state_id), "AdjointLayerThickness")
-                call get_option(trim(adj_sfield2%option_path) // "/prognostic/temporal_discretisation/theta", theta)
+                call get_option("/timestepping/theta", theta)
                 call get_option(trim(adj_sfield2%option_path) // "/prognostic/mean_layer_thickness", d0)
                 
                 big_mat => extract_block_csr_matrix(matrices, "InverseBigMatrix")

=== modified file 'assemble/Advection_Diffusion_DG.F90'
--- assemble/Advection_Diffusion_DG.F90	2012-10-11 17:52:01 +0000
+++ assemble/Advection_Diffusion_DG.F90	2012-10-23 15:04:25 +0000
@@ -1547,7 +1547,7 @@
              end if
           end if
        end if
-       
+
        ! Add stabilisation to the advection term if requested by the user.
        if (stabilisation_scheme==UPWIND) then
           ! NOTE: U_nl_q may (or may not) have been modified by the grid velocity
@@ -2603,7 +2603,6 @@
        nnAdvection_in=shape_shape(T_shape, T_shape_2, &
             &                       outer_advection_integral * detwei) 
 
-
     end if
     if (include_diffusion.or.have_buoyancy_adjustment_by_vertical_diffusion) then
        

=== added file 'assemble/Advection_Local_DG.F90'
--- assemble/Advection_Local_DG.F90	1970-01-01 00:00:00 +0000
+++ assemble/Advection_Local_DG.F90	2012-10-23 15:04:25 +0000
@@ -0,0 +1,2538 @@
+!    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 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 "fdebug.h"
+
+module advection_local_DG
+  !!< This module contains the Discontinuous Galerkin form of the advection
+  !!< equation when vector fields are represented in the local coordinates
+  !!< from the Piola transformation.
+  use elements
+  use global_parameters, only:current_debug_level, OPTION_PATH_LEN
+  use sparse_tools
+  use fetools
+  use dgtools
+  use fields
+  use fefields
+  use state_module
+  use shape_functions
+  use global_numbering
+  use transform_elements
+  use vector_tools
+  use fldebug
+  use vtk_interfaces
+  use Coordinates
+  use Tensors
+  use petsc_solve_state_module
+  use spud
+  use slope_limiters_dg
+  use sparsity_patterns
+  use sparse_matrices_fields
+  use sparsity_patterns_meshes
+  use Vector_Tools
+  use manifold_tools
+  use bubble_tools
+  use diagnostic_fields, only: calculate_diagnostic_variable
+  use global_parameters, only : FIELD_NAME_LEN
+
+  implicit none
+
+  ! Buffer for output messages.
+  character(len=255), private :: message
+
+  private
+  public solve_advection_dg_subcycle, solve_vector_advection_dg_subcycle,&
+       & solve_advection_cg_tracer, compute_U_residual, get_pv
+
+  !parameters specifying vector upwind options
+  integer, parameter :: VECTOR_UPWIND_EDGE=1, VECTOR_UPWIND_SPHERE=2
+
+  ! Local private control parameters. These are module-global parameters
+  ! because it would be expensive and/or inconvenient to re-evaluate them
+  ! on a per-element or per-face basis
+  real :: dt
+
+  ! Whether to include various terms
+  logical :: include_advection
+  contains
+
+    subroutine solve_advection_dg_subcycle(field_name, state, velocity_name,&
+         &continuity, flux)
+    !!< Construct and solve the advection equation for the given
+    !!< field using discontinuous elements.
+    
+    !! Name of the field to be solved for.
+    character(len=*), intent(in) :: field_name
+    !! Collection of fields defining system state.
+    type(state_type), intent(inout) :: state
+    !! Optional velocity name
+    character(len = *), intent(in) :: velocity_name
+    !! Solve continuity equation as opposed to advection equation
+    logical, intent(in), optional :: continuity
+    !! Return a flux variable such that T-T_old=div(Flux)
+    type(vector_field), intent(inout), optional :: Flux
+
+    !! Tracer to be solved for.
+    type(scalar_field), pointer :: T, T_old, s_field
+
+    !! Coordinate and velocity fields
+    type(vector_field), pointer :: X, U_nl
+
+    !! Velocity field in Cartesian coordinates for slope limiter
+    type(vector_field) :: U_nl_cartesian
+
+    !! Change in T over one timestep.
+    type(scalar_field) :: delta_T, delta_T_total
+
+    !! Sparsity of advection matrix.    
+    type(csr_sparsity), pointer :: sparsity
+    
+    !! System matrix.
+    type(csr_matrix) :: matrix, mass, inv_mass
+
+    !! Sparsity of mass matrix.
+    type(csr_sparsity) :: mass_sparsity
+
+    !! Right hand side vector.
+    type(scalar_field) :: rhs
+
+    !! Whether to invoke the slope limiter
+    logical :: limit_slope
+    !! Which limiter to use
+    integer :: limiter
+
+    !! Number of advection subcycles.
+    integer :: subcycles
+    real :: max_courant_number
+
+    !! Flux reconstruction stuff:
+    !! Upwind fluxes trace variable
+    type(scalar_field) :: UpwindFlux
+    !! Mesh where UpwindFlux lives
+    type(mesh_type), pointer :: lambda_mesh
+    !! Upwind Flux matrix
+    type(csr_matrix) :: UpwindFluxMatrix
+
+    character(len=FIELD_NAME_LEN) :: limiter_name
+    integer :: i, ele
+    real :: meancheck
+    ewrite(1,*) 'subroutine solve_advection_dg_subcycle'
+
+
+    T=>extract_scalar_field(state, field_name)
+    T_old=>extract_scalar_field(state, "Old"//field_name)
+    X=>extract_vector_field(state, "Coordinate")
+    U_nl=>extract_vector_field(state, velocity_name)
+
+    ewrite(2,*) 'UVALS in dg', maxval(abs(U_nl%val))
+
+    ! Reset T to value at the beginning of the timestep.
+    call set(T, T_old)
+    if(present(Flux)) then
+       call zero(Flux)
+    end if
+
+    sparsity => get_csr_sparsity_firstorder(state, T%mesh, T%mesh)
+    call allocate(matrix, sparsity) ! Add data space to the sparsity
+    ! pattern.
+    call zero(matrix)
+    
+    !Flux reconstruction allocations
+    if(present(Flux)) then
+       !Allocate trace variable for upwind fluxes
+       lambda_mesh=>extract_mesh(state, "VelocityMeshTrace")
+       call allocate(UpwindFlux,lambda_mesh, "UpwindFlux")
+       call zero(UpwindFlux)
+       !Allocate matrix that maps T values to upwind fluxes
+       sparsity => get_csr_sparsity_firstorder(state,lambda_mesh,T%mesh)
+       call allocate(UpwindFluxMatrix,sparsity)
+       call zero(UpwindFluxMatrix)
+    end if
+
+    mass_sparsity=make_sparsity_dg_mass(T%mesh)
+    call allocate(mass, mass_sparsity)
+    call zero(mass)
+    call allocate(inv_mass, mass_sparsity)
+    call zero(inv_mass)
+
+    ! Ensure delta_T inherits options from T.
+    call allocate(delta_T, T%mesh, "delta_T")
+    call zero(delta_T)
+    delta_T%option_path = T%option_path
+    if(present(Flux)) then
+       call allocate(delta_T_total, T%mesh, "deltaT_total")
+       call zero(delta_T_total)
+    end if
+    call allocate(rhs, T%mesh, trim(field_name)//" RHS")
+    call zero(rhs)
+
+    if(present(Flux)) then
+       call construct_advection_dg(matrix, rhs, field_name, state, &
+            mass, velocity_name=velocity_name,continuity=continuity, &
+            UpwindFlux=upwindflux, UpwindFluxMatrix=UpwindFluxMatrix)
+    else
+       call construct_advection_dg(matrix, rhs, field_name, state, &
+            mass, velocity_name=velocity_name,continuity=continuity)
+    end if
+
+    call get_dg_inverse_mass_matrix(inv_mass, mass)
+    
+    ! Note that since dt is module global, these lines have to
+    ! come after construct_advection_diffusion_dg.
+    call get_option("/timestepping/timestep", dt)
+    
+    if(have_option(trim(T%option_path)//&
+         &"/prognostic/temporal_discretisation/discontinuous_galerkin/&
+         &/number_advection_subcycles")) then
+       call get_option(trim(T%option_path)//&
+            &"/prognostic/temporal_discretisation/discontinuous_galerkin/&
+            &/number_advection_subcycles", subcycles)
+    else
+       call get_option(trim(T%option_path)//&
+            &"/prognostic/temporal_discretisation/discontinuous_galerkin/&
+            &maximum_courant_number_per_subcycle", Max_Courant_number)
+       
+       s_field => extract_scalar_field(state, "DG_CourantNumber")
+       call calculate_diagnostic_variable(state, "DG_CourantNumber_Local", &
+            & s_field)
+       ewrite_minmax(s_field)
+       subcycles = ceiling( maxval(s_field%val)/Max_Courant_number)
+       call allmax(subcycles)
+       ewrite(2,*) 'Number of subcycles for tracer eqn: ', subcycles
+    end if
+
+    limit_slope=.false.
+    if (have_option(trim(T%option_path)//"/prognostic/spatial_discretisation/&
+         &discontinuous_galerkin/slope_limiter")) then
+       limit_slope=.true.
+
+       call get_option(trim(T%option_path)//"/prognostic/spatial_discretisation/discontinuous_galerkin/slope_limiter/name",limiter_name)
+
+       select case(trim(limiter_name))
+       case("Vertex_Based")
+          limiter=LIMITER_VB
+       case default
+          FLAbort('No such limiter')
+       end select
+       
+       ! Note unsafe for mixed element meshes
+       if (element_degree(T,1)==0) then
+          FLExit("Slope limiters make no sense for degree 0 fields")
+       end if
+
+    end if
+
+    if (limit_slope) then
+       ! Filter wiggles from T
+       call limit_vb_manifold(state,T,delta_T)
+       if(present(flux)) then
+          call zero(upwindflux)
+          call mult(delta_T_total, mass, delta_T)
+          call update_flux(Flux, Delta_T_total, UpwindFlux, Mass)
+       end if
+   end if
+
+    do i=1, subcycles
+       ewrite(1,*) 'SUBCYCLE', i
+
+       ! dT = Advection * T
+       call mult(delta_T, matrix, T)
+       ! dT = dT + RHS
+       call addto(delta_T, RHS, -1.0)
+       if(present(Flux)) then
+          if(maxval(abs(RHS%val))>1.0e-8) then
+             FLAbort('Flux reconstruction doesn''t work if diffusion/bcs present at the moment')
+          end if
+       end if
+       call scale(delta_T, -dt/subcycles)
+       ! dT = M^(-1) dT
+       if(present(flux)) then
+          call mult(UpwindFlux, upwindfluxmatrix, T)
+          call scale(UpwindFlux,-dt/subcycles)
+          call update_flux(Flux, Delta_T, UpwindFlux, Mass)
+       end if
+       ! T = T + dt/s * dT
+       call dg_apply_mass(inv_mass, delta_T)
+       ewrite(2,*) 'Delta_T', maxval(abs(delta_T%val))
+       call addto(T, delta_T)
+       !Probably need to halo_update(Flux) as well       
+       call halo_update(T)
+
+       if (limit_slope) then
+          ! Filter wiggles from T
+          call limit_vb_manifold(state,T,delta_T)
+          if(present(flux)) then
+             call mult(delta_t_total, mass, delta_t)
+             call zero(UpwindFlux)
+             call update_flux(Flux, Delta_T_total, UpwindFlux, Mass)
+          end if
+       end if
+    end do
+
+    if(present(Flux)) then
+       if(have_option('/material_phase::Fluid/scalar_field::LayerThickness/p&
+            &rognostic/spatial_discretisation/debug')) then
+          do ele = 1, ele_count(Flux)
+             call check_flux(Flux,T,T_old,X,mass,ele)
+          end do
+       end if
+       call deallocate(UpwindFlux)
+       call deallocate(UpwindFluxMatrix)
+       call deallocate(delta_T_total)
+    end if
+
+    call deallocate(delta_T)
+    call deallocate(matrix)
+    call deallocate(mass)
+    call deallocate(inv_mass)
+    call deallocate(mass_sparsity)
+    call deallocate(rhs)
+
+  end subroutine solve_advection_dg_subcycle
+
+  subroutine check_flux(Flux,T,T_old,X,mass,ele)
+    type(vector_field), intent(in) :: Flux, X
+    type(scalar_field), intent(in) :: T,T_old
+    integer, intent(in) :: ele
+    type(csr_matrix), intent(in) :: mass
+    !
+    real, dimension(Flux%dim,ele_loc(Flux,ele)) :: Flux_vals
+    real, dimension(ele_ngi(T,ele)) :: Div_Flux_gi
+    real, dimension(ele_ngi(T,ele)) :: Delta_T_gi, T_gi
+    real, dimension(ele_loc(T,ele)) :: Delta_T_rhs, Div_Flux_rhs,T_rhs
+    integer :: dim1, loc
+    real, dimension(ele_ngi(X,ele)) :: detwei
+    real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J_mat
+    real :: residual
+    Delta_T_gi = ele_val_at_quad(T,ele)-ele_val_at_quad(T_old,ele)
+    T_gi = ele_val_at_quad(T,ele)
+    Flux_vals = ele_val(Flux,ele)
+
+    Div_Flux_gi = 0.
+    !dn is loc x ngi x dim    
+    do dim1 = 1, Flux%dim
+       do loc = 1, ele_loc(Flux,ele)
+          Div_Flux_gi = Div_Flux_gi + Flux%mesh%shape%dn(loc,:,dim1)*&
+               &Flux_vals(dim1,loc)
+       end do
+    end do
+
+    call compute_jacobian(ele_val(X,ele), ele_shape(X,ele),&
+         detwei=detwei,J=J_mat)
+    
+    Delta_T_rhs = shape_rhs(ele_shape(T,ele),Delta_T_gi*detwei)
+    T_rhs = shape_rhs(ele_shape(T,ele),T_gi*detwei)
+    Div_Flux_rhs = shape_rhs(ele_shape(T,ele),Div_Flux_gi*&
+         Flux%mesh%shape%quadrature%weight)
+    
+    residual = &
+         maxval(abs(Delta_T_rhs-Div_Flux_rhs))/max(1.0&
+         &,maxval(abs(T_rhs)))
+    if(residual>1.0e-6) then
+       ewrite(2,*) residual, maxval(abs(T_rhs))
+       FLExit('Flux residual error.')
+    end if
+
+  end subroutine check_flux
+
+  subroutine update_flux(Flux, Delta_T, UpwindFlux, Mass)
+    !! Subroutine to compute div-conforming Flux such that
+    !! div Flux = Delta_T
+    !! with Flux.n = UpwindFlux on element boundaries
+    !! and then to add to Flux
+    type(vector_field), intent(inout) :: Flux
+    type(scalar_field), intent(in) :: Delta_T
+    type(scalar_field), intent(in) :: UpwindFlux
+    type(csr_matrix), intent(in) :: Mass
+    !
+    integer :: ele, i
+    real :: Area
+    integer, dimension(:), pointer :: T_ele
+
+    !Checking the Flux is in local representation
+    assert(Flux%dim==mesh_dim(Flux))
+
+    do ele = 1, ele_count(Flux)
+       T_ele => ele_nodes(delta_T,ele)
+       area = 0.
+       do i = 1, size(T_ele)
+          area = area + &
+               sum(row_val(Mass,t_ele(i)))
+       end do
+       call update_flux_ele(Flux, Delta_T, UpwindFlux,ele,area)
+    end do
+
+  end subroutine update_flux
+
+  subroutine update_flux_ele(Flux, Delta_T, UpwindFlux,ele,area)
+    !! Subroutine to compute div-conforming Flux such that
+    !! div Flux = Delta_T
+    !! with Flux.n = UpwindFlux on element boundaries
+    type(vector_field), intent(inout) :: Flux
+    type(scalar_field), intent(in) :: Delta_T
+    type(scalar_field), intent(in) :: UpwindFlux
+    integer, intent(in) :: ele
+    real, intent(in) :: area
+    !
+    real, dimension(Flux%dim,ele_loc(Flux,ele)) :: Flux_vals
+    real, dimension(ele_loc(Delta_T,ele)) :: Delta_T_vals
+    real, dimension(Flux%dim*ele_loc(Flux,ele),&
+         Flux%dim*ele_loc(Flux,ele)) :: Flux_mat
+    real, dimension(Flux%dim*ele_loc(Flux,ele)) :: Flux_rhs
+    integer :: row, ni, ele2, face, face2,floc, i, dim1, dim2
+    integer, dimension(:), pointer :: neigh
+    type(constraints_type), pointer :: flux_constraint
+    type(element_type), pointer :: flux_shape, T_shape
+    real, dimension(Flux%dim,ele_loc(Delta_T,ele),&
+         &ele_loc(Flux,ele)) :: grad_mat
+    real, dimension(ele_loc(Delta_t,ele)) :: delta_t_val
+    real, dimension(ele_loc(flux,ele),ele_loc(flux,ele)) &
+         &:: flux_mass_mat
+    real, dimension(ele_loc(Delta_t,ele)) :: div_flux_val
+    real, dimension(ele_ngi(Delta_t,ele)) :: div_flux_gi
+    real :: total_flux, residual
+    !! Need to set up and solve equation for flux DOFs
+    !! Rows are as follows:
+    !! All of the normal fluxes through edges, numbered according to
+    !! UpwindFlux numbering
+    !! Inner product of Flux with grad basis equaling gradient of delta_T
+    !! plus boundary conditions
+    !! Inner product of solution with curl basis equals zero
+    !! inner product of solution with constraint basis equals zero
+
+    !! Debugging check:
+    !! \int \Delta T dV = \int div Flux dV = \int Flux.n dS
+    !!                                     = \int \hat{Flux}.\hat{n dS}
+    neigh => ele_neigh(Flux,ele)
+    total_flux = 0.
+    do ni = 1, size(neigh)
+       ele2 = neigh(ni)
+       face = ele_face(Flux,ele,ele2)
+       if(ele2>0) then
+          face2 = ele_face(Flux,ele2,ele)
+       else
+          face2 = face
+       end if
+       if(face>face2) then
+          total_flux=total_flux+sum(face_val(UpwindFlux,face))
+       else
+          total_flux=total_flux-sum(face_val(UpwindFlux,face))
+       end if
+    end do
+    residual = abs(total_flux-sum(ele_val(Delta_T,ele)))&
+         &/max(1.0,maxval(ele_val(Delta_T,ele)))/area
+    if(residual>1.0e-10) then
+       ewrite(0,*) 'Residual = ', residual
+       FLExit('Bad residual')
+    end if
+
+    Flux_rhs = 0.
+    Flux_mat = 0.
+    flux_shape => ele_shape(flux,ele)
+    T_shape => ele_shape(delta_t,ele)
+    flux_constraint => flux%mesh%shape%constraints
+
+    row = 1
+    !first do the face DOFs
+    neigh => ele_neigh(Flux,ele)
+    do ni = 1, size(neigh)
+       ele2 = neigh(ni)
+       face = ele_face(Flux,ele,ele2)
+       if(ele2>0) then
+          face2 = ele_face(Flux,ele2,ele)
+       else
+          face2 = face
+       end if
+       floc = face_loc(upwindflux,face)
+
+       call update_flux_face(&
+            Flux,UpwindFlux,face,face2,&
+            ele,ele2,Flux_mat(row:row+floc-1,:),&
+            Flux_rhs(row:row+floc-1))
+       row = row + floc
+    end do
+
+    !Next the gradient basis
+    !Start with \nabla\cdot F = \Delta T
+    !Integrate with test function over element
+    ! <\phi,Div F>_E = <\phi, \Delta T>_E
+
+    grad_mat = shape_dshape(T_shape, flux_shape%dn,&
+         & T_shape%quadrature%weight)
+    delta_t_val = ele_val(delta_t,ele)
+    do i = 1, ele_loc(Delta_T,ele)-1
+       do dim1 = 1, flux%dim
+          Flux_mat(row,(dim1-1)*ele_loc(flux,ele)+1:dim1*ele_loc(flux,ele))&
+               &=grad_mat(dim1,i,:)
+       end do
+       flux_rhs(row) = delta_t_val(i)
+       row = row + 1
+    end do
+
+    !Then the curl basis
+    !Adding in the curl terms
+    flux_mass_mat = shape_shape(flux_shape,flux_shape,T_shape%quadrature%weight)
+    do i = 1, flux_constraint%n_curl_basis
+       flux_mat(row,:) = 0.
+       flux_rhs(row) = 0.
+       do dim1 = 1, mesh_dim(flux)
+          do dim2 = 1, mesh_dim(flux)
+             flux_mat(&
+                  row,(dim2-1)*ele_loc(flux,ele)+1:dim2*ele_loc(flux,ele)) =&
+                  flux_mat(&
+                  row,(dim2-1)*ele_loc(flux,ele)+1:dim2*ele_loc(flux,ele)) +&
+                  matmul(flux_constraint%curl_basis(i,:,dim1),&
+                  flux_mass_mat)
+          end do
+       end do
+       row = row + 1
+    end do
+
+    !Then the constraints
+    do i = 1, flux_constraint%n_constraints
+       do dim1 = 1, mesh_dim(flux)
+          flux_mat(&
+               row,(dim1-1)*ele_loc(flux,ele)+1:&
+               dim1*ele_loc(flux,ele)) =&
+               flux_constraint%orthogonal(i,:,dim1)
+       end do
+       row = row + 1
+    end do
+
+    call solve(flux_mat,flux_rhs)
+
+    do dim1 = 1, mesh_dim(flux)
+       flux_vals(dim1,:) = &
+            flux_rhs((dim1-1)*ele_loc(flux,ele)+1:dim1*ele_loc(flux,ele))
+    end do
+
+    !A debugging check.
+    !Computing < \phi, div F> in local coordinates
+    !and comparing with <\phi, delta_T> in global coordinates
+    !(works because of pullback formula)
+
+    !dn is loc x ngi x dim
+    !Flux is dim x loc
+    div_flux_gi = 0.
+    do dim1 =1, size(flux_vals,1)
+       do i=1,size(flux_vals,2)
+          div_flux_gi = div_flux_gi + flux_vals(dim1,i)*&
+               &flux_shape%dn(i,:,dim1)
+       end do
+    end do
+    div_flux_val = shape_rhs(T_shape,div_flux_gi*T_shape%quadrature%weight)
+    residual = maxval(abs(div_flux_val-delta_T_val))/max(1.0&
+         &,maxval(abs(delta_T_val)))/area
+    assert(residual<1.0e-10)
+
+    !Add the flux in
+    call addto(flux,ele_nodes(flux,ele),flux_vals)
+  end subroutine update_flux_ele
+
+  subroutine update_flux_face(&
+       Flux,UpwindFlux,face,face2,&
+       ele,ele2,Flux_mat_rows,Flux_rhs_rows)
+    type(vector_field), intent(in) :: Flux
+    type(scalar_Field), intent(in) :: UpwindFlux
+    integer, intent(in) :: face,face2,ele,ele2
+    real, dimension(face_loc(upwindflux,ele),&
+         &mesh_dim(flux)*ele_loc(flux,ele)), &
+         &intent(inout) :: flux_mat_rows
+    real, dimension(flux%mesh%shape%constraints%n_face_basis), intent(inout)&
+         &:: flux_rhs_rows
+    !
+    integer :: dim1, row, flux_dim1
+    real, dimension(mesh_dim(flux), face_ngi(flux, face)) :: n_local
+    real, dimension(face_ngi(flux,face)) :: detwei_f
+    real :: weight
+    type(element_type), pointer :: flux_face_shape, upwindflux_face_shape
+    integer, dimension(face_loc(flux,face)) :: flux_face_nodes
+    real, dimension(mesh_dim(flux),face_loc(upwindflux,face),&
+         face_loc(flux,face)) :: face_mat
+    !
+    flux_face_shape=>face_shape(flux, face)
+    upwindflux_face_shape=>face_shape(upwindflux, face)
+    flux_face_nodes=face_local_nodes(flux, face)
+
+    !Get normal in local coordinates
+    call get_local_normal(n_local,weight,flux,&
+         &local_face_number(flux%mesh,face))
+
+    detwei_f = weight*flux_face_shape%quadrature%weight
+    face_mat = shape_shape_vector(&
+         upwindflux_face_shape,flux_face_shape,detwei_f,n_local)
+    !Equation is:
+    ! <\phi,Flux.n> = <\phi,UpwindFlux> for all trace test functions \phi.
+
+    do row = 1, size(flux_mat_rows,1)
+       flux_mat_rows(row,:) = 0.
+       do dim1 = 1,mesh_dim(flux)
+          flux_dim1 = (dim1-1)*ele_loc(flux,ele)
+          flux_mat_rows(row,flux_dim1+flux_face_nodes)&
+               &=face_mat(dim1,row,:)
+       end do
+    end do
+    !Need to adopt a sign convention:
+    !If face>face_2
+    !We store flux with normal pointing from face into face_2
+    !Otherwise 
+    !We store flux with normal pointing from face_2 into face
+    ! Outflow boundary integral.
+    if(face>face2) then
+       flux_rhs_rows = face_val(UpwindFlux,face)
+    else
+       flux_rhs_rows = -face_val(UpwindFlux,face)
+    end if
+  end subroutine update_flux_face
+
+  subroutine construct_advection_dg(big_m, rhs, field_name,&
+       & state, mass, velocity_name, continuity,&
+       & UpwindFlux, UpwindFluxMatrix)
+    !!< Construct the advection equation for discontinuous elements in
+    !!< acceleration form.
+    !!< 
+    !!< If mass is provided then the mass matrix is not added into big_m or
+    !!< rhs. It is instead returned as mass. This may be useful for testing
+    !!< or for solving equations otherwise than in acceleration form.
+
+    !! Main advection matrix.    
+    type(csr_matrix), intent(inout) :: big_m
+    !! Right hand side vector.
+    type(scalar_field), intent(inout) :: rhs
+
+    !! Name of the field to be advected.
+    character(len=*), intent(in) :: field_name
+    !! Collection of fields defining system state.
+    type(state_type), intent(inout) :: state
+    !! Optional separate mass matrix.
+    type(csr_matrix), intent(inout), optional :: mass
+    !! Optional velocity name
+    character(len = *), intent(in), optional :: velocity_name
+    !! solve the continuity equation
+    logical, intent(in), optional :: continuity
+    !! Upwind flux field
+    type(scalar_field), intent(in), optional :: UpwindFlux
+    !! Upwind flux matrix
+    type(csr_matrix), intent(inout), optional :: UpwindFluxMatrix
+
+    !! Position, and velocity fields.
+    type(vector_field) :: X, U, U_nl
+    !! Tracer to be solved for.
+    type(scalar_field) :: T
+
+    !! Local velocity name
+    character(len = FIELD_NAME_LEN) :: lvelocity_name
+
+    !! Element index
+    integer :: ele
+
+    !! Status variable for field extraction.
+    integer :: stat
+
+    ewrite(1,*) "Writing advection equation for "&
+         &//trim(field_name)
+
+    ! These names are based on the CGNS SIDS.
+    T=extract_scalar_field(state, field_name)
+    X=extract_vector_field(state, "Coordinate")
+
+    if(present(velocity_name)) then
+       lvelocity_name = velocity_name
+    else
+       lvelocity_name = "NonlinearVelocity"
+    end if
+
+    if (.not.have_option(trim(T%option_path)//"/prognostic&
+         &/spatial_discretisation/discontinuous_galerkin&
+         &/advection_scheme/none")) then
+       U_nl=extract_vector_field(state, lvelocity_name)
+       call incref(U_nl)
+       include_advection=.true.
+    else
+       ! Forcing a zero NonlinearVelocity will disable advection.
+       U=extract_vector_field(state, "Velocity", stat)
+       if (stat/=0) then 
+          FLExit("Oh dear, no velocity field. A velocity field is required for advection!")
+       end if
+       call allocate(U_nl, U%dim, U%mesh, "LocalNonlinearVelocity")
+       call zero(U_nl)
+       include_advection=.false.
+    end if
+
+    assert(has_faces(X%mesh))
+    assert(has_faces(T%mesh))
+
+    call zero(big_m)
+    call zero(RHS)
+    call zero(mass)
+
+    element_loop: do ele=1,element_count(T)
+
+       call construct_adv_element_dg(ele, big_m, rhs,&
+            & X, T, U_nl, mass, continuity=continuity,&
+            & upwindflux=upwindflux, upwindfluxmatrix=upwindfluxmatrix)
+
+    end do element_loop
+
+    ! Drop any extra field references.
+
+    call deallocate(U_nl)
+
+  end subroutine construct_advection_dg
+
+  subroutine construct_adv_element_dg(ele, big_m, rhs,&
+       & X, T, U_nl, mass, continuity, upwindflux, upwindfluxmatrix)
+    !!< Construct the advection_diffusion equation for discontinuous elements in
+    !!< acceleration form.
+    implicit none
+    !! Index of current element
+    integer :: ele
+    !! Main advection matrix.
+    type(csr_matrix), intent(inout) :: big_m
+    !! Right hand side vector.
+    type(scalar_field), intent(inout) :: rhs
+    !! Optional separate mass matrix.
+    type(csr_matrix), intent(inout) :: mass
+
+    !! Position and velocity.
+    type(vector_field), intent(in) :: X, U_nl
+
+    type(scalar_field), intent(in) :: T
+
+    !! Upwind flux field
+    type(scalar_field), intent(in), optional :: UpwindFlux
+    !! Upwind flux matrix
+    type(csr_matrix), intent(inout), optional :: UpwindFluxMatrix
+
+
+    !! Solve the continuity equation rather than advection
+    logical, intent(in), optional :: continuity
+
+    ! Bilinear forms.
+    real, dimension(ele_loc(T,ele), ele_loc(T,ele)) :: &
+         mass_mat
+    real, dimension(ele_loc(T,ele), ele_loc(T,ele)) :: &
+         Advection_mat, Ad_mat1, Ad_mat2
+
+    ! Local assembly matrices.
+    real, dimension(ele_loc(T,ele), ele_loc(T,ele)) :: l_T_mat
+    real, dimension(ele_loc(T,ele)) :: l_T_rhs
+
+    ! Local variables.
+
+    ! Neighbour element, face and neighbour face.
+    integer :: ele_2, face, face_2
+    ! Loops over faces.
+    integer :: ni
+
+    ! Transform from local to physical coordinates.
+    real, dimension(ele_ngi(X,ele)) :: detwei
+    real, dimension(mesh_dim(U_nl), X%dim, ele_ngi(X,ele)) :: J_mat 
+
+    ! Different velocities at quad points.
+    real, dimension(U_nl%dim, ele_ngi(U_nl, ele)) :: U_quad
+    real, dimension(U_nl%dim, ele_ngi(U_nl, ele)) :: U_nl_q
+    real, dimension(ele_ngi(U_nl, ele)) :: U_nl_div_q
+
+    ! Node and shape pointers.
+    integer, dimension(:), pointer :: T_ele
+    type(element_type), pointer :: T_shape, U_shape
+    ! Neighbours of this element.
+    integer, dimension(:), pointer :: neigh
+
+    integer :: gi,i,j
+
+    !----------------------------------------------------------------------
+    ! Establish local node lists
+    !----------------------------------------------------------------------
+
+    T_ele=>ele_nodes(T,ele)  ! Tracer node numbers
+
+    !----------------------------------------------------------------------
+    ! Establish local shape functions
+    !----------------------------------------------------------------------
+
+    T_shape=>ele_shape(T, ele)
+    U_shape=>ele_shape(U_nl, ele)
+
+    !==========================
+    ! Coordinates
+    !==========================
+
+    ! Get J_mat
+    call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), detwei=detwei, J=J_mat)
+
+    !----------------------------------------------------------------------
+    ! Construct bilinear forms.
+    !----------------------------------------------------------------------
+
+    ! Element density matrix.
+    !  /
+    !  | T T  dV
+    !  / 
+    mass_mat = shape_shape(T_shape, T_shape, detwei)
+
+    if (include_advection) then
+
+       ! Advecting velocity at quadrature points.
+       U_quad=ele_val_at_quad(U_nl,ele)
+       U_nl_q=U_quad
+
+       U_nl_div_q=ele_div_at_quad(U_nl, ele, U_shape%dn)
+
+       ! Element advection matrix
+       !    /                           /
+       !  - | (grad phi dot U_nl) T dV -| phi ( div U_nl ) T dV
+       !    /                           /
+       Advection_mat = -dshape_dot_vector_shape(T_shape%dn, U_nl_q, T_shape,&
+            & T_shape%quadrature%weight)
+       if(.not.present_and_true(continuity)) then
+          Advection_mat = Advection_mat&
+               &-shape_shape(T_shape, T_shape, U_nl_div_q * T_shape&
+               &%quadrature%weight)
+       end if
+    else
+       Advection_mat=0.0
+    end if
+
+    !----------------------------------------------------------------------
+    ! Perform global assembly.
+    !----------------------------------------------------------------------
+
+    l_T_rhs=0.0
+
+    ! Assemble matrix.
+
+    ! Advection.
+    l_T_mat= Advection_mat
+
+    call addto(mass, t_ele, t_ele, mass_mat)
+
+    call addto(big_m, t_ele, t_ele, l_T_mat)
+
+    !-------------------------------------------------------------------
+    ! Interface integrals
+    !-------------------------------------------------------------------
+
+    neigh=>ele_neigh(T, ele)
+
+    neighbourloop: do ni=1,size(neigh)
+
+       !----------------------------------------------------------------------
+       ! Find the relevant faces.
+       !----------------------------------------------------------------------
+
+       ! These finding routines are outside the inner loop so as to allow
+       ! for local stack variables of the right size in
+       ! construct_add_diff_interface_dg.
+
+       ele_2=neigh(ni)
+
+       ! Note that although face is calculated on field U, it is in fact
+       ! applicable to any field which shares the same mesh topology.
+       face=ele_face(T, ele, ele_2)
+
+       if (ele_2>0) then
+          ! Internal faces.
+          face_2=ele_face(T, ele_2, ele)
+       else
+          ! External face.
+          face_2=face
+       end if
+
+       call construct_adv_interface_dg(ele, face, face_2,&
+            & big_m, rhs, X, T, U_nl, upwindflux, upwindfluxmatrix)
+
+    end do neighbourloop
+
+  end subroutine construct_adv_element_dg
+
+  subroutine construct_adv_interface_dg(ele, face, face_2, &
+       big_m, rhs, X, T, U_nl,upwindflux, upwindfluxmatrix)
+
+    !!< Construct the DG element boundary integrals on the ni-th face of
+    !!< element ele.
+    implicit none
+
+    integer, intent(in) :: ele, face, face_2
+    type(csr_matrix), intent(inout) :: big_m
+    type(scalar_field), intent(inout) :: rhs
+    ! We pass these additional fields to save on state lookups.
+    type(vector_field), intent(in) :: X, U_nl
+    type(scalar_field), intent(in) :: T
+    !! Upwind flux field
+    type(scalar_field), intent(in), optional :: UpwindFlux
+    !! Upwind flux matrix
+    type(csr_matrix), intent(inout), optional :: UpwindFluxMatrix
+
+    ! Face objects and numberings.
+    type(element_type), pointer :: T_shape, T_shape_2
+    integer, dimension(face_loc(T,face)) :: T_face
+    integer, dimension(face_loc(T,face_2)) :: T_face_2
+    integer, dimension(face_loc(U_nl,face)) :: U_face
+    integer, dimension(face_loc(U_nl,face_2)) :: U_face_2
+
+    ! Note that both sides of the face can be assumed to have the same
+    ! number of quadrature points.
+    real, dimension(U_nl%dim, face_ngi(U_nl, face)) :: n1, n2, normal, U_nl_q,&
+         & u_f_q, u_f2_q, div_u_f_q
+    logical, dimension(face_ngi(U_nl, face)) :: inflow
+    real, dimension(face_ngi(U_nl, face)) :: U_nl_q_dotn1, U_nl_q_dotn2, U_nl_q_dotn, income
+    ! Variable transform times quadrature weights.
+    real, dimension(face_ngi(T,face)) :: detwei
+    real, dimension(U_nl%dim, X%dim, face_ngi(T,face)) :: J
+    real, dimension(face_ngi(T,face)) :: inner_advection_integral, outer_advection_integral
+
+    ! Bilinear forms
+    real, dimension(face_loc(T,face),face_loc(T,face)) :: nnAdvection_out
+    real, dimension(face_loc(T,face),face_loc(T,face_2)) :: nnAdvection_in
+
+    ! Normal weights
+    real :: w1, w2
+    integer :: i
+
+    T_face=face_global_nodes(T, face)
+    T_shape=>face_shape(T, face)
+
+    T_face_2=face_global_nodes(T, face_2)
+    T_shape_2=>face_shape(T, face_2)
+
+    !Get normals
+    call get_local_normal(n1, w1, U_nl, local_face_number(U_nl%mesh,face))
+    call get_local_normal(n2, w2, U_nl, local_face_number(U_nl%mesh,face_2))    
+
+    if (include_advection) then
+
+       ! Advecting velocity at quadrature points.
+       U_f_q = face_val_at_quad(U_nl, face)
+       U_f2_q = face_val_at_quad(U_nl, face_2)
+
+       u_nl_q_dotn1 = sum(U_f_q*w1*n1,1)
+       u_nl_q_dotn2 = -sum(U_f2_q*w2*n2,1)
+       U_nl_q_dotn=0.5*(u_nl_q_dotn1+u_nl_q_dotn2)
+
+       ! Inflow is true if the flow at this gauss point is directed
+       ! into this element.
+       inflow = u_nl_q_dotn<0.0
+       income = merge(1.0,0.0,inflow)
+
+       !----------------------------------------------------------------------
+       ! Construct bilinear forms.
+       !----------------------------------------------------------------------
+
+       ! Calculate outflow boundary integral.
+       ! can anyone think of a way of optimising this more to avoid
+       ! superfluous operations (i.e. multiplying things by 0 or 1)?
+
+       ! first the integral around the inside of the element
+       ! (this is the flux *out* of the element)
+       inner_advection_integral = (1.-income)*u_nl_q_dotn
+       nnAdvection_out=shape_shape(T_shape, T_shape,  &
+            &                     inner_advection_integral*T_shape%quadrature%weight)
+
+       ! now the integral around the outside of the element
+       ! (this is the flux *in* to the element)
+       outer_advection_integral = income*u_nl_q_dotn
+       nnAdvection_in=shape_shape(T_shape, T_shape_2, &
+            &                       outer_advection_integral*T_shape%quadrature%weight)
+
+    end if
+
+    !----------------------------------------------------------------------
+    ! Perform global assembly.
+    !----------------------------------------------------------------------
+
+    ! Insert advection in matrix.
+    if (include_advection) then
+
+       ! Outflow boundary integral.
+       call addto(big_M, T_face, T_face,&
+            nnAdvection_out)
+
+       ! Inflow boundary integral.
+       call addto(big_M, T_face, T_face_2,&
+            nnAdvection_in)
+
+       if(present(UpwindFlux).and.present(UpwindFluxMatrix)) then
+          call construct_upwindflux_interface(upwindflux,upwindfluxmatrix)
+       end if
+    end if
+
+  contains 
+    subroutine construct_upwindflux_interface(upwindflux,upwindfluxmatrix)
+      type(scalar_field), intent(in) :: upwindflux
+      type(csr_matrix), intent(inout) :: upwindfluxmatrix
+      !
+      ! Bilinear forms
+      real, dimension(face_loc(upwindflux,face),&
+           face_loc(T,face)) :: upwindflux_mat_in,upwindflux_mat_out
+      type(element_type), pointer :: flux_shape
+      integer, dimension(face_loc(upwindflux,face)) :: flux_face
+      integer :: sign
+
+      flux_shape=>face_shape(upwindflux, face)
+      flux_face=face_global_nodes(upwindflux, face)
+
+      upwindflux_mat_in=shape_shape(flux_shape, T_shape,  &
+           income*u_nl_q_dotn*T_shape%quadrature%weight)
+      upwindflux_mat_out=shape_shape(flux_shape, T_shape,  &
+           (1.0-income)*u_nl_q_dotn*T_shape%quadrature%weight)
+
+      !Need to adopt a sign convention:
+      !If face>face_2
+      !We store flux with normal pointing from face into face_2
+      !Otherwise 
+      !We store flux with normal pointing from face_2 into face
+      ! Outflow boundary integral.
+
+      !inflow = u_nl_q_dotn<0.0
+      !income = 1 if(inflow) and 0 otherwise
+      if(face>face_2) then
+          sign = 1
+       else
+          sign = -1
+       end if
+       
+       !Factor of 0.5 is because we visit from both sides
+       if(face.ne.face_2) then
+          call addto(upwindfluxmatrix, flux_face, T_face_2,&
+               0.5*sign*upwindflux_mat_in)
+       end if
+       call addto(upwindfluxmatrix, flux_face, T_face,&
+            0.5*sign*upwindflux_mat_out)
+
+      end subroutine construct_upwindflux_interface
+
+  end subroutine construct_adv_interface_dg  
+
+  subroutine solve_advection_cg_tracer(Q,D,D_old,Flux,U_nl,QF,state)
+    !!< Solve the continuity equation for D*Q with Flux
+    !!< d/dt (D*Q) + div(Flux*Q) = diffusion terms.
+    !!< Done on the vorticity mesh
+    !!< Return a PV flux Q defined on the velocity mesh
+    !!< which is the projection of Flux*Q into the velocity space
+    !!< Note that Flux and Q contain a factor of dt for convenience.
+    type(state_type), intent(inout) :: state
+    type(scalar_field), intent(in),target :: D,D_old
+    type(scalar_field), intent(inout),target :: Q
+    type(vector_field), intent(in) :: Flux,U_nl
+    type(vector_field), intent(inout) :: QF
+    !
+    type(csr_sparsity), pointer :: Q_sparsity
+    type(csr_matrix) :: Adv_mat
+    type(scalar_field) :: Q_rhs, Qtest1,Qtest2
+    type(vector_field), pointer :: X, down
+    type(scalar_field), pointer :: Q_old
+    integer :: ele
+    real :: dt, t_theta, residual
+
+    ewrite(1,*) '  subroutine solve_advection_cg_tracer('
+
+    X=>extract_vector_field(state, "Coordinate")
+    down=>extract_vector_field(state, "GravityDirection")
+    Q_old=>extract_scalar_field(state, "Old"//trim(Q%name))
+    call get_option('/timestepping/timestep',dt)
+    call get_option('/timestepping/theta',t_theta)
+    
+    !set up matrix and rhs
+    Q_sparsity => get_csr_sparsity_firstorder(state, Q%mesh, Q%mesh)
+    call allocate(adv_mat, Q_sparsity) ! Add data space to the sparsity
+    ! pattern.
+    call zero(adv_mat)
+    call allocate(Q_rhs,Q%mesh,trim(Q%name)//"RHS")
+    call zero(Q_rhs)
+    call zero(QF)
+
+    do ele = 1, ele_count(Q)
+       call construct_advection_cg_tracer_ele(Q_rhs,adv_mat,Q,D,D_old,Flux&
+            &,X,down,U_nl,dt,t_theta,ele)
+    end do
+
+    ewrite(2,*) 'Q_RHS', maxval(abs(Q_rhs%val))
+
+    call petsc_solve(Q,adv_mat,Q_rhs)
+
+    !! Compute the PV flux to pass to velocity equation
+    do ele = 1, ele_count(Q)
+       call construct_pv_flux_ele(QF,Q,Q_old,D,D_old,Flux,&
+            X,down,U_nl,t_theta,ele)
+    end do
+
+    if(have_option('/material_phase::Fluid/scalar_field::PotentialVorticity/&
+         &prognostic/debug')) then
+       call allocate(Qtest1,Q%mesh,trim(Q%name)//"test1")
+       call zero(Qtest1)
+       call allocate(Qtest2,Q%mesh,trim(Q%name)//"test2")
+       call zero(Qtest2)
+       do ele = 1, ele_count(Q)
+          call test_pv_flux_ele(Qtest1,Qtest2,QF,Q,Q_old,D,D_old,&
+               Flux,X,down,t_theta,ele)
+       end do
+       ewrite(2,*) 'Error = ', maxval(abs(Qtest2%val)), maxval(abs(Qtest1%val))
+       ewrite(2,*) 'Error = ', maxval(abs(Qtest1%val-Qtest2%val)), maxval(Qtest1%val)
+       residual = maxval(abs(Qtest1%val-Qtest2%val))/max(1.0&
+            &,maxval(abs(Qtest1%val)))
+       ewrite(2,*) 'residual from qtest', residual
+       assert(residual<1.0e-6)
+       ewrite(2,*) 'test passed'
+       call deallocate(Qtest1)
+       call deallocate(Qtest2)
+    end if
+
+    call deallocate(adv_mat)
+    call deallocate(Q_rhs)
+
+    ewrite(1,*) 'END  subroutine solve_advection_cg_tracer('
+
+  end subroutine solve_advection_cg_tracer
+  
+  subroutine construct_advection_cg_tracer_ele(Q_rhs,adv_mat,Q,D,D_old,Flux,&
+       & X,down,U_nl,dt,t_theta,ele)
+    type(scalar_field), intent(in) :: D,D_old,Q
+    type(scalar_field), intent(inout) :: Q_rhs
+    type(csr_matrix), intent(inout) :: Adv_mat
+    type(vector_field), intent(in) :: X, Flux, down,U_nl
+    integer, intent(in) :: ele
+    real, intent(in) :: dt, t_theta
+    !
+    real, dimension(ele_loc(Q,ele),ele_loc(Q,ele)) :: l_adv_mat,tmp_mat
+    real, dimension(ele_loc(Q,ele)) :: l_rhs
+    real, dimension(U_nl%dim,ele_ngi(U_nl,ele)) :: U_nl_gi
+    real, dimension(X%dim, ele_ngi(U_nl,ele)) :: U_cart_gi
+    real, dimension(Flux%dim,ele_ngi(Flux,ele)) :: Flux_gi
+    real, dimension(Flux%dim,FLux%dim,ele_ngi(Flux,ele)) :: Grad_Flux_gi,&
+         & tensor_gi
+    real, dimension(ele_ngi(X,ele)) :: detwei, Q_gi, D_gi, &
+         & D_old_gi,div_flux_gi, detwei_l, detJ
+    real, dimension(ele_loc(D,ele)) :: D_val
+    real, dimension(ele_loc(Q,ele)) :: Q_val
+    real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J
+    type(element_type), pointer :: Q_shape, Flux_shape
+    integer :: loc,dim1,dim2,gi
+    real :: tau, alpha, tol, area, h,c_sc,disc_indicator
+    real, dimension(ele_ngi(X,ele),ele_ngi(X,ele)) :: Metric, MetricT
+    real, dimension(mesh_dim(X),ele_ngi(X,ele)) :: gradQ
+    type(element_type), pointer :: D_shape
+    !
+
+    ! Get J and detwei
+    call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), detwei&
+         &=detwei, J=J, detJ=detJ)
+
+    D_shape => ele_shape(D,ele)
+    Q_shape => ele_shape(Q,ele)
+    Flux_shape => ele_shape(Flux,ele)
+    D_val = invert_pi_ele(ele_val(D,ele),D_shape,detwei)
+    D_gi = matmul(transpose(D_shape%n),D_val)
+    D_val = invert_pi_ele(ele_val(D_old,ele),D_shape,detwei)
+    D_old_gi = matmul(transpose(D_shape%n),D_val)
+    Q_gi = ele_val_at_quad(Q,ele)
+    Flux_gi = ele_val_at_quad(Flux,ele)
+    U_nl_gi = ele_val_at_quad(U_nl,ele)
+    Q_val = ele_val(Q,ele)
+
+    !Equations
+    ! D^{n+1} = D^n + \pi div Flux 
+    ! If define \int_K\phi\hat{D}/J J dS = \int_K\phi D J dS, then
+    ! <\phi,\pi(\hat{D}/J)> = <\phi, D> for all \phi,
+    ! => \pi \hat{D} = D
+    ! \hat{D}^{n+1}/J = \hat{D}^n/J + div Flux
+
+    ! So (integrating by parts)
+    ! <\gamma, \hat{D}^{n+1}/J> = 
+    !             <\gamma, D^n/J> - <\nabla\gamma, F>
+
+    ! and a consistent theta-method for Q is
+    ! <\gamma,\hat{D}^{n+1}Q^{n+1}/J> = 
+    !                  <\gamma, \hat{D}^nQ^n/J> - 
+    !                                     \theta<\nabla\gamma,FQ^{n+1}> 
+    !                                  -(1-\theta)<\nabla\gamma,FQ^n>
+
+    detwei_l = Q_shape%quadrature%weight
+    !Advection terms
+    !! <-grad gamma, FQ >
+    !! Can be evaluated locally using pullback
+    tmp_mat = dshape_dot_vector_shape(&
+         Q_shape%dn,Flux_gi,Q_shape,detwei_l)
+    l_adv_mat = t_theta*tmp_mat
+    l_rhs = -(1-t_theta)*matmul(tmp_mat,Q_val)
+
+    !Mass terms
+    !! <gamma, QD>
+    !! Requires transformed integral
+    tmp_mat = shape_shape(ele_shape(Q,ele),ele_shape(Q,ele),D_gi*detwei_l)
+    l_adv_mat = l_adv_mat + tmp_mat
+    tmp_mat = shape_shape(ele_shape(Q,ele),ele_shape(Q,ele),D_old_gi*detwei_l)
+    l_rhs = l_rhs + matmul(tmp_mat,Q_val)
+
+    !Laplacian filter options
+    if(have_option(trim(Q%option_path)//'/prognostic/spatial_discretisation/&
+         &continuous_galerkin/laplacian_filter')) then
+       FLExit('Laplacian filter not coded yet.')
+    end if
+    !Streamline upwinding options
+    if(have_option(trim(Q%option_path)//'/prognostic/spatial_discretisation/&
+         &continuous_galerkin/streamline_upwinding')) then
+       !! Replace test function \gamma with \gamma + \tau u\cdot\nabla\gamma
+       call get_option(trim(Q%option_path)//'/prognostic/spatial_discretisat&
+            &ion/continuous_galerkin/streamline_upwinding/alpha',alpha)
+       call get_option(trim(Q%option_path)//'/prognostic/spatial_discretisat&
+            &ion/continuous_galerkin/streamline_upwinding/tol',tol)
+
+       !Gradients of fluxes
+       Grad_flux_gi = ele_grad_at_quad(Flux,ele,Flux_shape%dn)
+       div_flux_gi = 0.
+       tensor_gi = 0.
+       do dim1 = 1, Flux%dim
+          div_flux_gi = div_flux_gi + grad_flux_gi(dim1,dim1,:)
+          do dim2 = 1, Flux%dim
+             tensor_gi(dim1,dim2,:) = U_nl_gi(dim1,:)*Flux_gi(dim2,:)
+          end do
+       end do
+       !real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J
+       do gi = 1, ele_ngi(X,ele)
+          U_cart_gi(:,gi) = matmul(transpose(J(:,:,gi)),U_nl_gi(:,gi))/detJ(gi)
+       end do
+
+       if(maxval(sqrt(sum(U_cart_gi**2,1)))<tol) then
+          tau = 0.
+       else
+          area = sum(detwei)
+          h = sqrt(4*area/sqrt(3.))
+          tau = h*alpha/maxval(sqrt(sum(U_cart_gi**2,1)))
+       end if
+       !! Mass terms
+       !! tau < u . grad gamma, QD>
+       !! = tau < grad gamma, u QD>
+       !! can do locally because of pullback formula
+       tmp_mat = dshape_dot_vector_shape(&
+            Q_shape%dn,U_nl_gi,Q_shape,D_gi*detwei_l/detJ)
+       l_adv_mat = l_adv_mat + tau*tmp_mat
+       tmp_mat = dshape_dot_vector_shape(&
+            Q_shape%dn,U_nl_gi,Q_shape,D_old_gi*detwei_l/detJ)
+       l_rhs = l_rhs + tau*matmul(tmp_mat,Q_val)
+       !! Advection terms
+       !! tau < u. grad gamma, div (F Q)>
+       !! = tau <grad gamma, u (F.grad Q + Q div F)>
+       !! can do locally because of pullback formula
+       !! tensor = u outer F (F already contains factor of dt)
+       tmp_mat = dshape_dot_vector_shape(&
+            Q_shape%dn,U_nl_gi,Q_shape,div_flux_gi*detwei_l/detJ) &
+            + dshape_tensor_dshape(&
+            Q_shape%dn,tensor_gi,Q_shape%dn,detwei_l/detJ)
+       l_adv_mat = l_adv_mat - tau*t_theta*tmp_mat
+       l_rhs = l_rhs + tau*(1-t_theta)*matmul(tmp_mat,Q_val)
+    end if
+    !Discontinuity capturing options
+    if(have_option(trim(Q%option_path)//'/prognostic/spatial_discretisation/&
+         &continuous_galerkin/discontinuity_capturing')) then
+       ! call have_option(&
+       !      &trim(Q%option_path)//'/prognostic/spatial_discretisation/&
+       !      &continuous_galerkin/discontinuity_capturing/&
+       !      &scaling_coefficient',c_sc)
+
+       ! do gi=1,ele_ngi(Q,ele)
+       !    Metric(:,:,gi)=matmul(J(:,:,gi), transpose(J(:,:,gi)))/detJ(gi)
+       ! end do
+       ! MetricT(1,1,:) = MetricT(2,2,:)
+       ! MetricT(2,2,:) = MetricT(1,1,:)
+       ! MetricT(1,2,:) =-MetricT(2,1,:)
+       ! MetricT(2,1,:) =-MetricT(1,2,:)
+
+       ! area = sum(detwei)
+       ! D_gi = ele_val_at_quad(D_old,ele)
+       ! tmp_mat = dshape_tensor_dshape(Q_shape%dn, &
+       !      area*D_gi*MetricT, Q_shape%dn,&
+       !      Q_shape%quadrature%weight)
+       ! l_rhs = l_rhs - dt*c_sc*matmul(tmp_mat,Q_val)
+       FLAbort('not ready yet')
+    end if
+
+    call addto(Q_rhs,ele_nodes(Q_rhs,ele),l_rhs)
+    call addto(adv_mat,ele_nodes(Q_rhs,ele),ele_nodes(Q_rhs,ele),&
+         l_adv_mat)    
+
+  end subroutine construct_advection_cg_tracer_ele
+
+  subroutine construct_pv_flux_ele(QFlux,Q,Q_old,D,D_old,&
+       Flux,X,down,U_nl,t_theta,ele)
+    type(scalar_field), intent(in) :: Q,Q_old,D,D_old
+    type(vector_field), intent(inout) :: QFlux
+    type(vector_field), intent(in) :: X, Flux, Down, U_nl
+    integer, intent(in) :: ele
+    real, intent(in) :: t_theta
+    !
+    real, dimension(mesh_dim(X),ele_loc(QFlux,ele)) :: QFlux_perp_rhs
+    real, dimension(Flux%dim,ele_ngi(Flux,ele)) :: Flux_gi, Flux_perp_gi,&
+         QFlux_gi
+    real, dimension(ele_ngi(X,ele)) :: Q_gi,Q_old_gi,D_gi,D_old_gi,Qbar_gi
+    real, dimension(ele_loc(D,ele)) :: D_val
+    real, dimension(X%dim, ele_ngi(Q, ele)) :: up_gi
+    real, dimension(mesh_dim(QFlux),mesh_dim(QFlux),ele_loc(QFlux,ele)&
+         &,ele_loc(QFlux,ele)) :: l_u_mat
+    real, dimension(mesh_dim(QFlux)*ele_loc(QFlux,ele),mesh_dim(QFlux)&
+         &*ele_loc(QFlux,ele)) :: solve_mat
+    real, dimension(ele_loc(Q,ele)) :: Q_test1_rhs, Q_test2_rhs
+    real, dimension(mesh_dim(Qflux)*ele_loc(QFlux,ele)) :: solve_rhs
+    real, dimension(mesh_dim(Qflux),ele_ngi(Qflux,ele)) :: &
+         contravariant_velocity_gi, velocity_perp_gi
+    type(element_type), pointer :: Q_shape, QFlux_shape, Flux_shape
+    integer :: loc, orientation,gi, dim1, dim2
+    real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J
+    real, dimension(ele_ngi(X,ele)) :: detwei, detJ, detwei_l, div_flux_gi
+    real, dimension(mesh_dim(Q),ele_ngi(Q,ele)) :: gradQ_gi
+    real, dimension(Flux%dim,FLux%dim,ele_ngi(Flux,ele)) :: Grad_Flux_gi
+    real, dimension(U_nl%dim,ele_ngi(U_nl,ele)) :: U_nl_gi
+    real, dimension(X%dim,ele_ngi(U_nl,ele)) :: U_cart_gi
+    real :: residual, alpha, tol, tau, area, h
+    type(element_type), pointer :: D_shape
+    !
+
+    D_shape => ele_shape(D,ele)
+
+    call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), J, detwei, detJ)
+
+    Q_shape => ele_shape(Q,ele)
+    Flux_shape => ele_shape(Flux,ele)
+    QFlux_shape => ele_shape(QFlux,ele)
+    Q_gi = ele_val_at_quad(Q,ele)
+    Q_old_gi = ele_val_at_quad(Q_old,ele)
+    Qbar_gi = t_theta*Q_gi + (1-t_theta)*Q_old_gi
+    D_val = invert_pi_ele(ele_val(D,ele),D_shape,detwei)
+    D_gi = matmul(transpose(D_shape%n),D_val)
+    D_val = invert_pi_ele(ele_val(D_old,ele),D_shape,detwei)
+    D_old_gi = matmul(transpose(D_shape%n),D_val)
+    Flux_gi = ele_val_at_quad(Flux,ele)
+    U_nl_gi = ele_val_at_quad(U_nl,ele)
+
+    detwei_l = Q_shape%quadrature%weight
+    up_gi = -ele_val_at_quad(down,ele)
+
+    call get_up_gi(X,ele,up_gi,orientation)
+
+    !Equations
+    ! <\nabla \gamma, (-v,u)> = <(\gamma_x, \gamma_y),(-v,u)>
+    !                         = <(\gamma_y,-\gamma_x),( u,v)>
+    !                         = <-\nabla^\perp\gamma,(u,v)>
+
+    ! and a consistent theta-method for Q is
+    ! <\gamma, D^{n+1}Q^{n+1}> = <\gamma, D^nQ^n> - 
+    !                                     \theta<\nabla\gamma,FQ^{n+1}> 
+    !                                  -(1-\theta)<\nabla\gamma,FQ^n>
+    ! So vorticity update is 
+    ! <\gamma, \zeta^{n+1}> = <-\nabla^\perp\gamma,u^{n+1}>
+    != <-\nabla^\perp\gamma,u^n> 
+    !  +\theta<-\nabla^\perp\gamma,(FQ^{n+1})^\perp>
+    !  +(1-\theta)<-\nabla^\perp\gamma,(FQ^n)^\perp>
+    ! taking w = -\nabla^\perp\gamma,
+    ! <w,u^{n+1}>=<w,u^n> + <w,F(\theta Q^{n+1}+(1-\theta)Q^n)^\perp>
+    ! <w,u^{n+1}>=<w,u^n> + <w,\bar{Q}^\perp>
+    ! We actually store the inner product with test function (FQ_rhs)
+    ! (avoids having to do a solve)
+    
+    do gi  = 1, ele_ngi(QFlux,ele)
+       QFlux_gi(:,gi) = Flux_gi(:,gi)*Qbar_gi(gi)
+    end do
+
+    !! Additional dissipative terms.
+    !! They all take the form
+    !! <\nabla gamma, G> 
+    !! Which can be evaluated in local coordinates due to pullback magic
+
+    !Laplacian filter options
+    if(have_option(trim(Q%option_path)//'/prognostic/spatial_discretisation/&
+         &continuous_galerkin/laplacian_filter')) then
+       FLExit('Laplacian filter not coded yet.')
+    end if
+    !Streamline upwinding options
+    if(have_option(trim(Q%option_path)//'/prognostic/spatial_discretisation/&
+         &continuous_galerkin/streamline_upwinding')) then
+       !! Replace test function \gamma with \gamma + \tau u\cdot\nabla\gamma
+       call get_option(trim(Q%option_path)//'/prognostic/spatial_discretisat&
+            &ion/continuous_galerkin/streamline_upwinding/alpha',alpha)
+       call get_option(trim(Q%option_path)//'/prognostic/spatial_discretisat&
+            &ion/continuous_galerkin/streamline_upwinding/tol',tol)
+
+       !Gradients of fluxes
+       Grad_flux_gi = ele_grad_at_quad(Flux,ele,Flux_shape%dn)
+       div_flux_gi = 0.
+       do dim1 = 1, Flux%dim
+          div_flux_gi = div_flux_gi + grad_flux_gi(dim1,dim1,:)
+       end do
+       !real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J
+       do gi = 1, ele_ngi(X,ele)
+          U_cart_gi(:,gi) = matmul(transpose(J(:,:,gi)),U_nl_gi(:,gi))/detJ(gi)
+       end do
+
+       if(maxval(sqrt(sum(U_cart_gi**2,1)))<tol) then
+          tau = 0.
+       else
+          area = sum(detwei)
+          h = sqrt(4*area/sqrt(3.))
+          tau = h*alpha/maxval(sqrt(sum(U_cart_gi**2,1)))
+       end if
+       !! Mass terms
+       !! tau < u . grad gamma, -(QD)^{n+1}+(QD)^n>
+       !! = tau < grad gamma, u(-(QD)^{n+1}+(QD)^n)>
+       !! can do locally because of pullback formula
+       !! Extra flux is tau u\Delta(QD)
+       !! minus sign from definition of vorticity
+       !! tau
+       do gi = 1, ele_ngi(Q,ele)
+          QFLux_gi(:,gi) = QFlux_gi(:,gi) + &
+               & tau*U_nl_gi(:,gi)*(Q_gi(gi)*D_gi(gi)-Q_old_gi(gi)&
+               &*D_old_gi(gi))/detJ(gi)
+
+       end do
+       !! Advection terms
+       !! tau < u. grad gamma, div (F\bar{Q})>
+       !! can do locally because of pullback formula
+       !! = <grad gamma, tau u (F.grad\bar{Q} + \bar{Q} div F)/det J>_{Ehat}
+       !! minus sign because of definition of vorticity
+
+       GradQ_gi = t_theta*ele_grad_at_quad(Q,ele,Q_shape%dn) + &
+            (1-t_theta)*ele_grad_at_quad(Q_old,ele,Q_shape%dn)
+
+       do gi = 1, ele_ngi(Q,ele)
+          QFlux_gi(:,gi) = QFlux_gi(:,gi) - tau*u_nl_gi(:,gi)*(&
+               & sum(Flux_gi(:,gi)*gradQ_gi(:,gi)) + &
+               & div_flux_gi(gi)*Qbar_gi(gi))/detJ(gi)
+       end do
+    end if
+    !Discontinuity capturing options
+    if(have_option(trim(Q%option_path)//'/prognostic/spatial_discretisation/&
+         &continuous_galerkin/discontinuity_capturing')) then
+       FLExit('Discontinuity capturing not coded yet.')
+    end if
+
+    !! Evaluate 
+    !! < w, F^\perp > in local coordinates
+
+    Flux_perp_gi(1,:) = -orientation*QFlux_gi(2,:)
+    Flux_perp_gi(2,:) =  orientation*QFlux_gi(1,:)
+
+    QFlux_perp_rhs = shape_vector_rhs(QFlux_shape,Flux_perp_gi,&
+         & QFlux_shape%quadrature%weight)
+
+    call set(QFlux,ele_nodes(QFlux,ele),QFlux_perp_rhs)
+
+  end subroutine construct_pv_flux_ele
+
+  subroutine test_pv_flux_ele(Qtest1,Qtest2,QFlux,Q,Q_old,D,D_old,&
+       Flux,X,down,t_theta,ele)
+    type(scalar_field), intent(in) :: Q,Q_old,D,D_old
+    type(scalar_field), intent(inout) :: Qtest1, Qtest2
+    type(vector_field), intent(in) :: X, Flux, Down,QFlux
+    integer, intent(in) :: ele
+    real, intent(in) :: t_theta
+    !
+    real, dimension(mesh_dim(X),ele_loc(QFlux,ele)) :: QFlux_perp_rhs
+    real, dimension(ele_ngi(X,ele)) :: Q_gi,Q_old_gi,D_gi,D_old_gi
+    real, dimension(Flux%dim,ele_ngi(Flux,ele)) :: Flux_gi, Flux_perp_gi
+    real, dimension(X%dim, ele_ngi(Q, ele)) :: up_gi
+    real, dimension(mesh_dim(QFlux), mesh_dim(QFlux), ele_ngi(QFlux,ele)) &
+         :: Metric, Metricf
+    real, dimension(mesh_dim(QFlux),mesh_dim(QFlux),ele_loc(QFlux,ele)&
+         &,ele_loc(QFlux,ele)) :: l_u_mat
+    real, dimension(mesh_dim(QFlux)*ele_loc(QFlux,ele),mesh_dim(QFlux)&
+         &*ele_loc(QFlux,ele)) :: solve_mat
+    real, dimension(ele_loc(Q,ele)) :: Q_test1_rhs, Q_test2_rhs
+    real, dimension(mesh_dim(Qflux)*ele_loc(QFlux,ele)) :: solve_rhs
+    real, dimension(mesh_dim(Qflux),ele_ngi(Qflux,ele)) :: &
+         contravariant_velocity_gi, velocity_perp_gi
+    type(element_type), pointer :: Q_shape, QFlux_shape,D_shape
+    integer :: loc, orientation,gi, dim1, dim2
+    real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J
+    real, dimension(ele_ngi(X,ele)) :: detwei, detJ
+    real, dimension(X%dim, X%dim, ele_ngi(X,ele)) :: rot
+    real, dimension(ele_loc(D,ele)) :: D_val
+
+    D_shape => ele_shape(D,ele)
+    Q_shape => ele_shape(Q,ele)
+    QFlux_shape => ele_shape(QFlux,ele)
+    Q_gi = ele_val_at_quad(Q,ele)
+    Q_old_gi = ele_val_at_quad(Q_old,ele)
+
+    call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), J, detwei, detJ)
+    D_val = invert_pi_ele(ele_val(D,ele),D_shape,detwei)
+    D_gi = matmul(transpose(D_shape%n),D_val)
+    D_val = invert_pi_ele(ele_val(D_old,ele),D_shape,detwei)
+    D_old_gi = matmul(transpose(D_shape%n),D_val)
+
+    Flux_gi = ele_val_at_quad(Flux,ele)
+    Qflux_perp_rhs = ele_val(Qflux,ele)
+
+    do gi  = 1, ele_ngi(QFlux,ele)
+       Flux_gi(:,gi) = Flux_gi(:,gi)*(t_theta*Q_gi(gi) + &
+         (1-t_theta)*Q_old_gi(gi))
+    end do
+
+    up_gi = -ele_val_at_quad(down,ele)
+    call get_up_gi(X,ele,up_gi,orientation)
+
+    do gi=1, ele_ngi(QFlux,ele)
+       rot(1,:,gi)=(/0.,-up_gi(3,gi),up_gi(2,gi)/)
+       rot(2,:,gi)=(/up_gi(3,gi),0.,-up_gi(1,gi)/)
+       rot(3,:,gi)=(/-up_gi(2,gi),up_gi(1,gi),0./)
+    end do
+    do gi=1,ele_ngi(QFlux,ele)
+       Metric(:,:,gi)=matmul(J(:,:,gi), transpose(J(:,:,gi)))/detJ(gi)
+       Metricf(:,:,gi)=matmul(J(:,:,gi), &
+            matmul(rot(:,:,gi), transpose(J(:,:,gi))))/detJ(gi)
+    end do
+
+    l_u_mat = shape_shape_tensor(qflux_shape, qflux_shape, &
+         qflux_shape%quadrature%weight,Metric)
+    solve_mat = 0.
+    solve_rhs = 0.
+    do dim1 = 1, mesh_dim(qflux)
+       do dim2 = 1, mesh_dim(qflux)
+          solve_mat( (dim1-1)*ele_loc(qflux,ele)+1:dim1*ele_loc(qflux,ele),&
+               (dim2-1)*ele_loc(qflux,ele)+1:dim2*ele_loc(qflux,ele)) = &
+               l_u_mat(dim1,dim2,:,:)
+       end do
+       solve_rhs( (dim1-1)*ele_loc(qflux,ele)+1:dim1*ele_loc(qflux,ele))&
+            = QFlux_perp_rhs(dim1,:)
+    end do
+    call solve(solve_mat,solve_rhs)
+    !Don't need to include constraints since u eqn is
+    ! <w,Q^\perp> + <<[w],\lambda>> + \sum_i C_i\cdot w \Gamma_i = RHS
+    ! but if w=-\nabla^\perp\gamma then [w]=0 and C_i\cdot w=0.
+
+    do dim1 = 1, mesh_dim(qflux)
+       flux_perp_gi(dim1,:) = matmul(transpose(QFlux_shape%n),&
+            solve_rhs((dim1-1)*ele_loc(qflux,ele)+1:dim1*ele_loc(qflux&
+            &,ele)))
+    end do
+    !Now compute vorticity
+    do gi=1,ele_ngi(X,ele)
+       Metric(:,:,gi)=matmul(J(:,:,gi), transpose(J(:,:,gi)))/detJ(gi)
+       contravariant_velocity_gi(:,gi) = &
+            matmul(Metric(:,:,gi),Flux_perp_gi(:,gi))
+    end do
+
+    !!Q_test1_rhs contains -<\nabla\gamma,Q>
+
+    select case(mesh_dim(X))
+    case (2)
+       velocity_perp_gi(1,:) = -contravariant_velocity_gi(2,:)
+       velocity_perp_gi(2,:) =  contravariant_velocity_gi(1,:)
+       Q_test1_rhs = dshape_dot_vector_rhs(Q%mesh%shape%dn, &
+            velocity_perp_gi,Q%mesh%shape%quadrature%weight)
+       Q_test1_rhs = q_test1_rhs*orientation
+    case default
+       FLAbort('Exterior derivative not implemented for given mesh dimension')
+    end select
+
+    !!Q_test2_rhs contains <\gamma,Q^{n+1}\tilde{D}^{n+1}-
+    !!                             Q^n\tilde{D}^n>
+
+    Q_test2_rhs = shape_rhs(ele_shape(Q,ele),(Q_gi*D_gi-Q_old_gi&
+         &*D_old_gi)*Q_shape%quadrature%weight)
+
+    call addto(Qtest1,ele_nodes(Q,ele),Q_test1_rhs)
+    call addto(Qtest2,ele_nodes(Q,ele),Q_test2_rhs)
+
+  end subroutine test_pv_flux_ele
+
+  subroutine get_PV(state,PV,velocity,D,path)
+    type(state_type), intent(inout) :: state
+    type(vector_field), intent(in) :: velocity
+    type(scalar_field), intent(inout) :: PV
+    type(scalar_field), intent(in) :: D
+    character(len=OPTION_PATH_LEN), optional :: path
+    !
+    type(vector_field), pointer :: X, down
+    type(scalar_field), pointer :: Coriolis
+    type(csr_matrix) :: pv_mass_matrix
+    type(csr_sparsity), pointer :: pv_mass_sparsity, curl_sparsity
+    type(scalar_field) :: pv_rhs
+    integer :: ele
+    !!DEbugging
+    integer ::dim1
+    real, dimension(3,3) :: X_val
+
+    ewrite(1,*) 'subroutine get_pv'
+
+    pv_mass_sparsity => &
+         get_csr_sparsity_firstorder(state, pv%mesh, pv%mesh)
+    call allocate(pv_mass_matrix,pv_mass_sparsity)
+    call zero(pv_mass_matrix)
+
+    X=>extract_vector_field(state, "Coordinate")
+    down=>extract_vector_field(state, "GravityDirection")
+    Coriolis=>extract_scalar_field(state, "Coriolis")
+
+    call allocate(pv_rhs, pv%mesh, 'PvRHS')
+    call zero(pv_rhs)
+    
+    do ele = 1, ele_count(pv)
+       call assemble_pv_ele(pv_mass_matrix,pv_rhs,velocity,D,Coriolis,X&
+            &,down,ele)
+    end do
+
+    if(present(path)) then
+       call petsc_solve(pv,pv_mass_matrix,pv_rhs,option_path=path)
+    else
+       call petsc_solve(pv,pv_mass_matrix,pv_rhs)
+    end if
+
+    call deallocate(pv_mass_matrix)
+    call deallocate(pv_rhs)
+
+  end subroutine get_pv
+
+  subroutine assemble_PV_ele(pv_mass_matrix,pv_rhs,velocity,D,Coriolis,X&
+       &,down,ele)
+    !!Subroutine to compute the pv from a *local velocity* field
+    type(vector_field), intent(in) :: X, down, velocity
+    type(scalar_field), intent(inout) :: pv_rhs
+    type(scalar_field), intent(in) :: D, Coriolis
+    type(csr_matrix), intent(inout) :: pv_mass_matrix
+    integer, intent(in) :: ele
+    !
+    real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J
+    real, dimension(ele_ngi(X,ele)) :: detwei, detJ, l_f, l_d
+    real, dimension(ele_loc(D,ele)) :: D_val
+    real, dimension(ele_loc(pv_rhs,ele),ele_loc(pv_rhs,ele))&
+         :: l_mass_mat
+    real, dimension(ele_loc(pv_rhs,ele))&
+         :: l_rhs
+    real, dimension(mesh_dim(velocity),ele_ngi(velocity,ele)) :: velocity_gi,&
+         contravariant_velocity_gi, velocity_perp_gi
+    real, dimension(X%dim, ele_ngi(X,ele)) :: up_gi
+    integer :: orientation, gi
+    real, dimension(mesh_dim(X), mesh_dim(X), ele_ngi(X,ele)) :: Metric
+    type(element_type), pointer :: D_shape
+    !
+
+    D_shape => ele_shape(D,ele)
+
+    up_gi = -ele_val_at_quad(down,ele)
+    call get_up_gi(X,ele,up_gi,orientation)
+
+    call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), J, detwei, detJ)    
+    D_val = ele_val(D,ele)
+    D_val = invert_pi_ele(D_val,D_shape,detwei)
+    l_d = matmul(transpose(D_shape%n),D_val)
+
+    l_mass_mat = shape_shape(ele_shape(pv_rhs,ele),&
+         ele_shape(pv_rhs,ele),l_d*d_shape%quadrature%weight)
+
+    velocity_gi = ele_val_at_quad(velocity,ele)
+    do gi=1,ele_ngi(X,ele)
+       Metric(:,:,gi)=matmul(J(:,:,gi), transpose(J(:,:,gi)))/detJ(gi)
+       contravariant_velocity_gi(:,gi) = &
+            matmul(Metric(:,:,gi),velocity_gi(:,gi))
+    end do
+
+    !Relative vorticity
+
+    ! pv = \nabla^\perp\cdot u = v_x - u_y
+    ! < \gamma, pv > = <\gamma, v_x - u_y> = <-\gamma_x,v>+<\gamma_y,u>
+     ! < -\nabla^\perp \gamma, J^TJ u/detJ> in local coordinates
+    ! < \nabla \gamma, (J^TJ u)^\perp/detJ> in local coordinates
+    ! requires us to know the orientation of the manifold
+
+    select case(mesh_dim(X))
+    case (2)
+       velocity_perp_gi(1,:) = -contravariant_velocity_gi(2,:)
+       velocity_perp_gi(2,:) =  contravariant_velocity_gi(1,:)
+       l_rhs = dshape_dot_vector_rhs(pv_rhs%mesh%shape%dn, &
+            velocity_perp_gi,X%mesh%shape%quadrature%weight)
+       l_rhs = l_rhs*orientation
+    case default
+       FLAbort('Exterior derivative not implemented for given mesh dimension')
+    end select
+
+    ! Coriolis term
+    l_f = ele_val_at_quad(Coriolis,ele)
+    l_rhs = l_rhs + shape_rhs(ele_shape(pv_rhs,ele),l_f*detwei)
+
+    call addto(pv_mass_matrix,ele_nodes(pv_rhs,ele),&
+         ele_nodes(pv_rhs,ele),l_mass_mat)
+    call addto(pv_rhs,ele_nodes(pv_rhs,ele),l_rhs)
+
+  end subroutine assemble_pv_ele
+
+  function invert_pi_ele(D_val_in,D_shape,detwei) result (D_val_out)
+    real, intent(in), dimension(:) :: D_val_in
+    type(element_type), intent(in) :: D_shape
+    real, intent(in), dimension(:) :: detwei
+    real, dimension(size(D_val_in)) :: D_val_out !The result
+    !
+    real, dimension(size(D_val_in),size(D_val_in)) :: proj_mat, &
+         & mass_mat
+
+    !! \sum_i \phi_i \hat{D}_i w_i = \sum_i \phi_i D_i J_i w_i.
+
+    proj_mat = shape_shape(D_shape,D_shape,D_shape%quadrature%weight)
+    mass_mat = shape_shape(D_shape,D_shape,detwei)
+    D_val_out = matmul(mass_mat,D_val_in)
+    call solve(proj_mat,D_val_out)
+
+  end function invert_pi_ele
+
+  subroutine compute_U_residual(UResidual,oldU,oldD,newU,newD,PVFlux,state)
+    !!< Compute the residual in the U equation, given PVFlux computed
+    !!< from the PV advection equation. Note that the PVFlux has 
+    !!< already been perped, and multiplied by test functions.
+    !!< Residual will be multiplied by test functions
+    type(vector_field), intent(inout) :: UResidual
+    type(vector_field), intent(in) :: oldU,newU,PVFlux
+    type(scalar_field), intent(in) :: oldD,newD
+    type(state_type), intent(inout) :: state
+    !
+    integer :: ele, stat
+    type(vector_field), pointer :: X
+    type(scalar_field), pointer :: Orography
+    real :: dt, theta,g 
+    type(scalar_field), target :: l_orography
+    logical :: have_orography
+
+    ewrite(1,*) '  subroutine compute_U_residual('
+
+    X=>extract_vector_field(state, "Coordinate")
+    call get_option("/timestepping/timestep", dt)
+    call get_option("/timestepping/theta",theta)
+    call get_option("/physical_parameters/gravity/magnitude", g)
+
+    Orography=>extract_scalar_field(state, "Orography", stat)
+    have_orography = (stat==0)
+    if(.not.have_orography) then
+       call allocate(l_orography,X%mesh,"L_Orography")
+       orography => l_orography
+       call zero(l_orography)
+    end if
+
+    call zero(UResidual)
+    do ele = 1, ele_count(UResidual)
+       call compute_U_residual_ele(UResidual,oldU,oldD,newU,newD,PVFlux,&
+            orography,X,theta&
+            &,dt,g,ele)
+    end do
+    
+    if(.not.have_orography) then
+       call deallocate(l_orography)
+    end if
+
+    ewrite(1,*) 'END subroutine compute_U_residual('
+
+  end subroutine compute_U_residual
+
+  subroutine compute_U_residual_ele(UResidual,oldU,oldD,newU,newD,PVFlux,&
+       orography,X,theta&
+       &,dt,g,ele)
+    type(vector_field), intent(inout) :: UResidual
+    type(vector_field), intent(in) :: oldU,newU,PVFlux,X
+    type(scalar_field), intent(in) :: oldD,newD,orography
+    real, intent(in) :: dt,theta,g
+    integer, intent(in) :: ele
+    !
+    real, dimension(ele_ngi(X,ele)) :: detwei, detJ
+    real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J
+    real, dimension(mesh_dim(oldU), mesh_dim(oldU)) :: Metric
+    real, dimension(mesh_dim(oldU), ele_ngi(X,ele)) :: U_rhs, newU_rhs
+    real, dimension(ele_ngi(X,ele)) :: D_gi, newD_gi, D_bar_gi, K_bar_gi
+    real, dimension(mesh_dim(X), ele_ngi(X,ele)) :: U_gi,newU_gi
+    real, dimension(X%dim, ele_ngi(X,ele)) :: U_cart_gi,newU_cart_gi
+    real, dimension(mesh_dim(oldU),ele_loc(UResidual,ele)) :: UR_rhs
+    integer :: gi
+    type(element_type), pointer :: U_shape
+    real, dimension(ele_ngi(X,ele)) :: orography_gi
+
+    !r[w] = <w,u^{n+1}-u^n> - <w,FQ^\perp> 
+    !           - dt*<div w, g\bar{h} + \bar{|u|^2/2}>
+
+    call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), &
+         detwei=detwei, J=J, detJ=detJ)
+    U_shape=>ele_shape(oldU, ele)
+
+    !Get all the variables at the quadrature points
+    D_gi = ele_val_at_quad(oldD,ele)
+    orography_gi = ele_val_at_quad(orography,ele)
+    newD_gi = ele_val_at_quad(newD,ele)
+    U_gi = ele_val_at_quad(oldU,ele)
+    newU_gi = ele_val_at_quad(newU,ele)
+    U_rhs = 0.
+    newU_rhs = 0.
+    do gi=1,ele_ngi(oldU,ele)
+       Metric=matmul(J(:,:,gi), transpose(J(:,:,gi)))/detJ(gi)
+       U_rhs(:,gi) = matmul(Metric,U_gi(:,gi))
+       newU_rhs(:,gi) = matmul(Metric,newU_gi(:,gi))
+    end do
+    U_cart_gi = 0.
+    newU_cart_gi = 0.
+    do gi = 1, ele_ngi(oldU,ele)
+       U_cart_gi(:,gi) = matmul(transpose(J(:,:,gi)),U_gi(:,gi))/detJ(gi)
+       newU_cart_gi(:,gi) = matmul(transpose(J(:,:,gi)),newU_gi(:,gi))/detJ(gi)
+    end do
+
+    !!Residual is
+    !! U_new - (U_old + dt*Q^\perp - dt grad(g\bar{D} + \bar{K}))
+    !First the PV Flux (perped, contains factor of dt)
+    UR_rhs = -ele_val(PVFlux,ele)
+    !Now the time derivative term
+    UR_rhs = UR_rhs + shape_vector_rhs(U_shape,newU_rhs-U_rhs,&
+         U_shape%quadrature%weight)
+    !Now the gradient terms (done in local coordinates)
+    D_bar_gi = theta*newD_gi + (1-theta)*D_gi + orography_gi
+    K_bar_gi = 0.5*(theta*sum(newU_cart_gi**2,1)+(1-theta)*sum(U_cart_gi**2,1))
+    !integration by parts, so minus sign
+    UR_rhs = UR_rhs - dt * dshape_rhs(U_shape%dn,&
+         (g*D_bar_gi + K_bar_gi)*U_shape%quadrature%weight)
+
+    call set(UResidual,ele_nodes(UResidual,ele),UR_rhs)
+
+  end subroutine compute_U_residual_ele
+
+  subroutine solve_vector_advection_dg_subcycle(field_name, state, velocity_name)
+    !!< Construct and solve the advection equation for the given
+    !!< field using discontinuous elements.
+    
+    !! Name of the field to be solved for.
+    character(len=*), intent(in) :: field_name
+    !! Collection of fields defining system state.
+    type(state_type), intent(inout) :: state
+    !! Optional velocity name
+    character(len = *), intent(in) :: velocity_name
+
+    !! Velocity to be solved for.
+    type(vector_field), pointer :: U, U_old, U_cartesian
+
+    !! Coordinate and advecting velocity fields
+    type(vector_field), pointer :: X, U_nl, down
+
+    !! Temporary velocity fields
+    type(vector_field) :: U_tmp, U_cartesian_tmp
+
+    !! Velocity components Cartesian coordinates for slope limiter
+    type(scalar_field) :: U_component
+
+    !! Change in U over one timestep.
+    type(vector_field) :: delta_U, delta_U_tmp
+
+    !! DG Courant number field
+    type(scalar_field), pointer :: s_field
+
+    !! Sparsity of advection matrix.    
+    type(csr_sparsity), pointer :: sparsity
+    
+    !! System matrices.
+    type(block_csr_matrix) :: A, L, mass_local, inv_mass_local, &
+         inv_mass_cartesian
+
+    !! Sparsity of mass matrix.
+    type(csr_sparsity) :: mass_sparsity
+
+    !! Right hand side vector.
+    type(vector_field) :: rhs
+
+    !! Whether to invoke the slope limiter
+    logical :: limit_slope
+    !! Which limiter to use
+    integer :: limiter
+
+    !! Number of advection subcycles.
+    integer :: subcycles
+    real :: max_courant_number
+
+    character(len=FIELD_NAME_LEN) :: limiter_name
+    integer :: i, j, dim, ele, dim1
+    integer :: upwinding_option
+
+    if(have_option('/material_phase::Fluid/vector_field::Velocity/prognostic&
+         &/spatial_discretisation/discontinuous_galerkin/advection_scheme/ed&
+         &ge_coordinates_upwind')) then
+       upwinding_option = VECTOR_UPWIND_EDGE
+    else
+       if(have_option('/material_phase::Fluid/vector_field::Velocity/prognostic&
+            &/spatial_discretisation/discontinuous_galerkin/advection_scheme/sphere_coordinates_upwind')) then
+          upwinding_option = VECTOR_UPWIND_SPHERE
+       else
+          FLAbort('Unknown upwinding option')
+       end if
+    end if
+
+    U=>extract_vector_field(state, field_name)
+    U_old=>extract_vector_field(state, "Old"//field_name)
+    X=>extract_vector_field(state, "Coordinate")
+    U_cartesian=>extract_vector_field(state, "Velocity")
+    down=>extract_vector_field(state, "GravityDirection")
+
+    dim=mesh_dim(U)
+
+    ! Reset U to value at the beginning of the timestep.
+    call set(U, U_old)
+
+    sparsity => get_csr_sparsity_firstorder(state, U%mesh, U%mesh)
+
+    ! Add data space to the sparsity pattern.
+    call allocate(A, sparsity, (/dim,X%dim/))
+    call zero(A)
+    call allocate(L, sparsity, (/dim,X%dim/))
+    call zero(L)
+
+    mass_sparsity=make_sparsity_dg_mass(U%mesh)
+    call allocate(mass_local, mass_sparsity, (/dim,dim/))
+    call zero(mass_local)
+    call allocate(inv_mass_local, mass_sparsity, (/dim,dim/))
+    call zero(inv_mass_local)
+    call allocate(inv_mass_cartesian, mass_sparsity, (/X%dim,X%dim/))
+    call zero(inv_mass_cartesian)
+
+    ! Ensure delta_U inherits options from U.
+    call allocate(delta_U, U%dim, U%mesh, "delta_U")
+    call zero(delta_U)
+    call allocate(delta_U_tmp, U%dim, U%mesh, "delta_U")
+    call zero(delta_U_tmp)
+    delta_U%option_path = U_cartesian%option_path
+    call allocate(rhs, U%dim, U%mesh, trim(field_name)//" RHS")
+    call zero(rhs)
+
+    ! allocate temp velocity fields
+    call allocate(U_tmp, U%dim, U%mesh, 'tmpU')
+    call zero(U_tmp)
+    call allocate(U_cartesian_tmp, U_cartesian%dim, U_cartesian%mesh, 'tmpUcartesian')
+    call zero(U_cartesian_tmp)
+   
+    call construct_vector_advection_dg(A, L, mass_local, &
+         inv_mass_local, inv_mass_cartesian,&
+         rhs, field_name, state, upwinding_option, down, &
+         velocity_name=velocity_name)
+    
+    ! Note that since dt is module global, these lines have to
+    ! come after construct_advection_diffusion_dg.
+    call get_option("/timestepping/timestep", dt)
+    
+    if(have_option(trim(U_cartesian%option_path)//&
+         &"/prognostic/temporal_discretisation/discontinuous_galerkin/&
+         &/number_advection_subcycles")) then
+       call get_option(trim(U_cartesian%option_path)//&
+            &"/prognostic/temporal_discretisation/discontinuous_galerkin/&
+            &/number_advection_subcycles", subcycles)
+    else
+       call get_option(trim(U_cartesian%option_path)//&
+            &"/prognostic/temporal_discretisation/discontinuous_galerkin/&
+            &/maximum_courant_number_per_subcycle", Max_Courant_number)
+       
+       s_field => extract_scalar_field(state, "DG_CourantNumber")
+       call calculate_diagnostic_variable(state, "DG_CourantNumber_Local", &
+            & s_field)
+       ewrite_minmax(s_field)
+       
+       subcycles = ceiling( maxval(s_field%val)/Max_Courant_number)
+       call allmax(subcycles)
+       ewrite(2,*) 'Number of subcycles for velocity eqn: ', subcycles
+    end if
+
+    limit_slope=.false.
+    if (have_option(trim(U_cartesian%option_path)//&
+         "/prognostic/spatial_discretisation/&
+         &discontinuous_galerkin/slope_limiter")) then
+       limit_slope=.true.
+       
+       ! Note unsafe for mixed element meshes
+       if (element_degree(U,1)==0) then
+          FLExit("Slope limiters make no sense for degree 0 fields")
+       end if
+
+    end if
+
+    do i=1, subcycles
+       call mult_t(U_cartesian_tmp, L, U)
+       call mult(U_cartesian, inv_mass_cartesian, U_cartesian_tmp)
+        ! dU = Advection * U
+       ! A maps from cartesian to local
+       call mult(delta_U_tmp, A, U_cartesian)
+       ! dU = dU + RHS
+       !------------------------------------------
+       !Is there anything in RHS?
+       call addto(delta_U_tmp, RHS, -1.0)
+       !------------------------------------------
+       ! dU = M^(-1) dU
+       call mult(delta_U, inv_mass_local, delta_U_tmp)
+       !------------------------------------------
+       !Do we need the below?
+       call mult_t(U_cartesian_tmp, L, delta_U)
+       call mult(U_cartesian, inv_mass_cartesian, U_cartesian_tmp)
+       !------------------------------------------
+
+       ! U = U + dt/s * dU
+       call addto(U, delta_U, scale=-dt/subcycles)
+       call halo_update(U)
+
+       if (limit_slope) then
+
+          call mult_t(U_cartesian_tmp, L, U)
+          call mult(U_cartesian, inv_mass_cartesian, U_cartesian_tmp)
+
+          !Really crap limiter
+          !Just limit each cartesian component
+          do dim1 = 1, U_cartesian%dim
+             u_component = extract_scalar_field(U_cartesian, dim1)
+             call limit_vb(state, U_component)
+          end do
+          !limiter=VECTOR_LIMITER_VB
+          !call limit_slope_dg(U_cartesian, state, limiter)
+
+          call mult(U_tmp, L, U_cartesian)
+          call mult(U, inv_mass_local, U_tmp)
+       end if
+
+    end do
+
+    call deallocate(delta_U)
+    call deallocate(A)
+    call deallocate(L)
+    call deallocate(mass_local)
+    call deallocate(inv_mass_local)
+    call deallocate(inv_mass_cartesian)
+    call deallocate(mass_sparsity)
+    call deallocate(rhs)
+
+  end subroutine solve_vector_advection_dg_subcycle
+
+  subroutine construct_vector_advection_dg(A, L, mass_local, inv_mass_local,&
+       & inv_mass_cartesian, rhs, field_name, &
+       & state, upwinding_option, down, velocity_name) 
+    !!< Construct the advection equation for discontinuous elements in
+    !!< acceleration form.
+    !!< 
+    !!< If mass is provided then the mass matrix is not added into big_m or
+    !!< rhs. It is instead returned as mass. This may be useful for testing
+    !!< or for solving equations otherwise than in acceleration form.
+
+    !! Main advection matrix.    
+    type(block_csr_matrix), intent(inout) :: A, L
+    !! Mass matrices.
+    type(block_csr_matrix), intent(inout) :: mass_local, inv_mass_local, &
+         inv_mass_cartesian
+    !! Right hand side vector.
+    type(vector_field), intent(inout) :: rhs
+    type(vector_field), intent(in) :: down
+    !! Name of the field to be advected.
+    character(len=*), intent(in) :: field_name
+    !! Collection of fields defining system state.
+    type(state_type), intent(inout) :: state
+    integer, intent(in) :: upwinding_option
+
+    !! Optional velocity name
+    character(len = *), intent(in), optional :: velocity_name
+
+    !! Position, velocity and advecting velocity fields.
+    type(vector_field) :: X, U, U_nl
+
+    !! Local velocity name
+    character(len = FIELD_NAME_LEN) :: lvelocity_name
+
+    !! Element index
+    integer :: ele
+
+    !! Status variable for field extraction.
+    integer :: stat
+
+    ewrite(1,*) "Writing advection equation for "&
+         &//trim(field_name)
+
+    ! These names are based on the CGNS SIDS.
+    U=extract_vector_field(state, field_name)
+    X=extract_vector_field(state, "Coordinate")
+
+    if(present(velocity_name)) then
+      lvelocity_name = velocity_name
+    else
+      lvelocity_name = "NonlinearVelocity"
+    end if
+
+    U_nl=extract_vector_field(state, lvelocity_name)
+    call incref(U_nl)
+
+    assert(has_faces(X%mesh))
+    assert(has_faces(U%mesh))
+    
+    call zero(A)
+    call zero(RHS)
+    call zero(mass_local)
+
+    element_loop: do ele=1,element_count(U)
+       
+       call construct_vector_adv_element_dg(ele, A, L,&
+            & mass_local, inv_mass_local, inv_mass_cartesian, rhs,&
+            & X, U, U_nl, upwinding_option, down)
+       
+    end do element_loop
+    
+    ! Drop any extra field references.
+
+    call deallocate(U_nl)
+
+  end subroutine construct_vector_advection_dg
+
+  subroutine construct_vector_adv_element_dg(ele, A, L,&
+       & mass_local, inv_mass_local, inv_mass_cartesian, rhs,&
+       & X, U, U_nl, upwinding_option, down)
+    !!< Construct the advection_diffusion equation for discontinuous elements in
+    !!< acceleration form.
+    implicit none
+    !! Index of current element
+    integer :: ele
+    !! Main advection matrix.
+    type(block_csr_matrix), intent(inout) :: A
+    !! Transformation matrix
+    type(block_csr_matrix), intent(inout) :: L
+    !! Mass matrices.
+    type(block_csr_matrix), intent(inout) :: mass_local, inv_mass_local, &
+         inv_mass_cartesian
+    !! Right hand side vector.
+    type(vector_field), intent(inout) :: rhs
+    
+    !! Position, velocity and advecting velocity.
+    type(vector_field), intent(in) :: X, U, U_nl, down
+
+    !! Upwinding option
+    integer, intent(in) :: upwinding_option
+
+    ! Bilinear forms.
+    real, dimension(mesh_dim(U),mesh_dim(U),ele_loc(U,ele),ele_loc(U,ele)) :: local_mass_mat
+    real, dimension(mesh_dim(U),X%dim,ele_loc(U,ele),ele_loc(U,ele)) :: l_mat
+    real, dimension(ele_loc(U,ele),ele_loc(U,ele)) :: cartesian_mass_mat
+    real, dimension(mesh_dim(U),X%dim,ele_loc(U,ele),ele_loc(U,ele)) :: &
+         Advection_mat
+
+    ! Local assembly matrices.
+    real, dimension(mesh_dim(U)*ele_loc(U,ele), mesh_dim(U)*ele_loc(U,ele)) :: l_mass, l_mass_cartesian
+
+    ! Local variables.
+    
+    ! Neighbour element, face and neighbour face.
+    integer :: ele_2, face, face_2
+    ! Loops over faces.
+    integer :: ni
+    
+    ! Transform from local to physical coordinates.
+    real, dimension(ele_ngi(X,ele)) :: detwei, detJ
+    real, dimension(mesh_dim(U), X%dim, ele_ngi(X,ele)) :: J, J_scaled
+    real, dimension(X%dim, mesh_dim(U), ele_ngi(X,ele)) :: pinvJ
+    real, dimension(ele_loc(U,ele), ele_ngi(X,ele), X%dim) :: dshape
+    real, dimension(mesh_dim(U), mesh_dim(U), ele_ngi(X,ele)) :: G
+
+    ! Different velocities at quad points.
+    real, dimension(U_nl%dim, ele_ngi(U_nl, ele)) :: U_nl_q
+    real, dimension(X%dim, ele_ngi(U_nl, ele)) :: U_cartesian_q
+    real, dimension(ele_ngi(U_nl, ele)) :: U_nl_div_q
+
+    ! Node and shape pointers.
+    integer, dimension(:), pointer :: U_ele
+    type(element_type), pointer :: U_shape, U_nl_shape
+    ! Neighbours of this element.
+    integer, dimension(:), pointer :: neigh
+
+    integer :: i, gi, dim, dim1, dim2, nloc, k
+
+    dim=mesh_dim(U)
+
+    !----------------------------------------------------------------------
+    ! Establish local node lists
+    !----------------------------------------------------------------------
+    
+    U_ele=>ele_nodes(U,ele)  ! Velocity node numbers
+
+    !----------------------------------------------------------------------
+    ! Establish local shape functions
+    !----------------------------------------------------------------------
+
+    U_shape=>ele_shape(U, ele)
+    U_nl_shape=>ele_shape(U_nl, ele)
+
+    !==========================
+    ! Coordinates
+    !==========================
+
+    ! Get J and detJ
+    call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), detwei=detwei, J=J, detJ=detJ)
+
+    !----------------------------------------------------------------------
+    ! Construct bilinear forms.
+    !----------------------------------------------------------------------
+
+    ! Element mass matrix.
+    !  /
+    !  | W G U  dV
+    !  / 
+    do gi=1,ele_ngi(X,ele)
+       G(:,:,gi)=matmul(J(:,:,gi),transpose(J(:,:,gi)))/detJ(gi)
+       J_scaled(:,:,gi)=J(:,:,gi)/detJ(gi)
+    end do
+
+    local_mass_mat=shape_shape_tensor(u_shape, u_shape, &
+         u_shape%quadrature%weight, G)
+
+    cartesian_mass_mat=shape_shape(u_shape, u_shape, detwei)
+
+    ! Transformation matrix
+    l_mat=shape_shape_tensor(u_shape, u_shape, u_shape%quadrature%weight, J)
+
+    ! Advecting velocity at quadrature points.
+    U_nl_q=ele_val_at_quad(U_nl,ele)
+    U_nl_div_q=ele_div_at_quad(U_nl, ele, U_shape%dn)
+
+    ! WE HAVE INTEGRATED BY PARTS TWICE
+    ! Element advection matrix
+    !    / 
+    !    | w . div (\bar{U} \tensor u )dV
+    !    / 
+    
+    ! becomes
+    
+    !    /
+    !    | \phi_i (dx/d\xi)^T.\div_L (\bar{U}_L \tensor \phi_j)dV_L
+    !    /
+    ! underscore U indicates local coordinates
+
+    ! split into two terms
+
+    !    /
+    !    | \phi_i (dx/d\xi)^T(\div_L \bar{U}_L)\phi_j dV_L
+    !    /
+
+    ! and 
+
+    !    /
+    !    | \phi_i (dx/d\xi)^T\bar{U}_L . grad_L \phi_j dV_L
+    !    /
+
+
+    U_cartesian_q=0.
+    do gi=1,ele_ngi(X,ele)
+       U_cartesian_q(:,gi)=matmul(transpose(J(:,:,gi)),U_nl_q(:,gi))
+       U_cartesian_q(:,gi)=U_cartesian_q(:,gi)/detJ(gi)
+    end do
+    Advection_mat = shape_vector_dot_dshape_tensor(U_shape, U_nl_q, U_shape&
+         &%dn, J_scaled, U_shape%quadrature%weight)
+
+   !----------------------------------------------------------------------
+   ! Perform global assembly.
+   !----------------------------------------------------------------------
+
+   ! Assemble matrices.
+
+   ! Return local mass separately.
+   do dim1 = 1,dim
+      do dim2 = 1,dim
+         call addto(mass_local, dim1, dim2, U_ele, U_ele, &
+              local_mass_mat(dim1,dim2,:,:))
+      end do
+   end do
+
+   ! calculate inverse_mass_local
+   nloc=ele_loc(U,ele)
+   do dim1 = 1,dim
+      do dim2 = 1,dim
+         l_mass(nloc*(dim1-1)+1:nloc*dim1, nloc*(dim2-1)+1:nloc*dim2)=&
+              local_mass_mat(dim1,dim2,:,:)
+      end do
+   end do
+
+   call invert(l_mass)
+
+   do dim1 = 1,dim
+      do dim2 = 1,dim
+         call addto(inv_mass_local, dim1, dim2, U_ele, U_ele, &
+              l_mass(nloc*(dim1-1)+1:nloc*dim1,nloc*(dim2-1)+1:nloc*dim2))
+      end do
+   end do
+
+   ! calculate inverse_mass_cartesian
+   call invert(cartesian_mass_mat)
+
+   do dim1 = 1,X%dim
+         call addto(inv_mass_cartesian, dim1, dim1, U_ele, U_ele, cartesian_mass_mat)
+   end do
+
+   ! advection matrix
+   do dim1 = 1,dim
+      do dim2 = 1,X%dim
+         call addto(A, dim1, dim2, U_ele, U_ele, Advection_mat(dim1,dim2,:,:))
+      end do
+   end do
+
+   ! transformation matrix
+   do dim1 = 1,dim
+      do dim2 = 1,X%dim
+         call addto(L, dim1, dim2, U_ele, U_ele, l_mat(dim1,dim2,:,:))
+      end do
+   end do
+
+   !-------------------------------------------------------------------
+   ! Interface integrals
+   !-------------------------------------------------------------------
+    
+   neigh=>ele_neigh(U, ele)
+
+   neighbourloop: do ni=1,size(neigh)
+
+      !----------------------------------------------------------------------
+      ! Find the relevant faces.
+      !----------------------------------------------------------------------
+       
+      ! These finding routines are outside the inner loop so as to allow
+      ! for local stack variables of the right size in
+      ! construct_add_diff_interface_dg.
+
+      ele_2=neigh(ni)
+       
+      ! Note that although face is calculated on field U, it is in fact
+      ! applicable to any field which shares the same mesh topology.
+      face=ele_face(U, ele, ele_2)
+    
+      if (ele_2>0) then
+         ! Internal faces.
+         face_2=ele_face(U, ele_2, ele)
+      else
+         ! External face.
+         face_2=face
+      end if
+
+      call construct_vector_adv_interface_dg(ele, ele_2, face, face_2,&
+           & A, rhs, X, U, U_nl, upwinding_option, down)
+
+   end do neighbourloop
+    
+ end subroutine construct_vector_adv_element_dg
+  
+  subroutine construct_vector_adv_interface_dg(ele, ele_2, face, face_2, &
+       & A, rhs, X, U, U_nl, upwinding_option, down)
+
+    !!< Construct the DG element boundary integrals on the ni-th face of
+    !!< element ele.
+    implicit none
+
+    integer, intent(in) :: ele, ele_2, face, face_2
+    type(block_csr_matrix), intent(inout) :: A
+    type(vector_field), intent(inout) :: rhs
+    ! We pass these additional fields to save on state lookups.
+    type(vector_field), intent(in) :: X, U, U_nl, down
+    integer, intent(in) :: upwinding_option
+
+    ! Face objects and numberings.
+    type(element_type), pointer :: U_shape, U_shape_2
+    integer, dimension(face_loc(U,face)) :: U_face, U_face_l
+    integer, dimension(face_loc(U,face_2)) :: U_face_2
+    integer, dimension(face_loc(U_nl,face)) :: U_nl_face
+    integer, dimension(face_loc(U_nl,face_2)) :: U_nl_face_2
+
+    ! Note that both sides of the face can be assumed to have the same
+    ! number of quadrature points.
+    real, dimension(U_nl%dim, face_ngi(U_nl, face)) :: n1_l, n2_l, U_nl_q,&
+         & U_nl_f_q, U_nl_f2_q
+    logical, dimension(face_ngi(U_nl, face)) :: inflow
+    real, dimension(face_ngi(U_nl, face)) :: U_nl_q_dotn1_l, U_nl_q_dotn2_l, U_nl_q_dotn_l, income
+    real, dimension(X%dim, face_ngi(X,face)) :: n1, n2, t11, t12, t21,t22,&
+         & pole_axis
+    real, dimension(X%dim, ele_ngi(X,ele)) :: up_gi, up_gi2
+    real, dimension(X%dim, face_ngI(X,face)) :: X_quad
+    real, dimension(X%dim, X%dim, face_ngi(X,face)) :: Btmp
+    real, dimension(mesh_dim(U), X%dim, face_ngi(X,face)) :: B
+    ! Variable transform times quadrature weights.
+    real, dimension(ele_ngi(X,ele)) :: detwei, detJ
+    real, dimension(face_ngi(X,face)) :: detwei_f, detJ_f
+    real, dimension(mesh_dim(U), X%dim, ele_ngi(X,ele)) :: J
+    real, dimension(mesh_dim(U), X%dim, face_ngi(X,face)) :: J_scaled
+    real, dimension(face_ngi(X,face)) :: inner_advection_integral, outer_advection_integral
+
+    ! Bilinear forms
+    real, dimension(mesh_dim(U),X%dim,face_loc(U,face),face_loc(U,face)) :: nnAdvection_out
+    real, dimension(mesh_dim(U),X%dim,face_loc(U,face),face_loc(U,face_2)) :: nnAdvection_in
+
+    ! normal weights
+    real :: w1, w2
+    integer :: dim, dim1, dim2, gi, i, k
+
+    dim=mesh_dim(U)
+
+    U_face=face_global_nodes(U, face)
+    U_shape=>face_shape(U, face)
+    X_quad = face_val_at_quad(X, face)
+
+    U_face_2=face_global_nodes(U, face_2)
+    U_shape_2=>face_shape(U, face_2)
+    
+    !Unambiguously calculate the normal using the face with the higher
+    !face number. This is so that the normal is identical on both sides.
+
+    call get_local_normal(n1_l, w1, U, local_face_number(U%mesh,face))
+    call get_local_normal(n2_l, w2, U, local_face_number(U%mesh,face_2))
+    
+    !----------------------------------------------------------------------
+    ! Construct element-wise quantities.
+    !----------------------------------------------------------------------
+
+    ! Advecting velocity at quadrature points.
+    U_nl_f_q = face_val_at_quad(U_nl, face)
+    U_nl_f2_q = face_val_at_quad(U_nl, face_2)
+
+    u_nl_q_dotn1_l = sum(U_nl_f_q*w1*n1_l,1)
+    u_nl_q_dotn2_l = -sum(U_nl_f2_q*w2*n2_l,1)
+    U_nl_q_dotn_l=0.5*(u_nl_q_dotn1_l+u_nl_q_dotn2_l)
+
+    ! Inflow is true if the flow at this gauss point is directed
+    ! into this element.
+    inflow = u_nl_q_dotn_l<0.0
+    income = merge(1.0,0.0,inflow)
+
+    ! Calculate tensor on face
+    if(ele_loc(X,ele).ne.3) then
+       ewrite(1,*) 'Nloc=',ele_loc(X,ele)
+       FLAbort('Hard coded for linear elements at the moment')
+    end if
+    call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), detwei=detwei, J=J, detJ=detJ)
+    forall(gi=1:face_ngi(X,face))
+       J_scaled(:,:,gi)=J(:,:,1)/detJ(1)
+    end forall
+
+    select case(upwinding_option)
+    case (VECTOR_UPWIND_EDGE)
+       ! Get outward pointing normals in physical space for face1
+       ! inward pointing normals in physical space for face2
+       n1=get_face_normal_manifold(X, ele, face)
+       !needs checking
+       if(ele_2<0) then
+          ! external boundary
+          n2=n1
+       else
+          n2=-get_face_normal_manifold(X, ele_2, face_2)
+       end if
+       ! Form 'bending' tensor
+       ! This rotates the normal to face 2 into -normal from face 1
+       ! u_b = t(t.u_b) + n_b(n_b.u_b)
+       ! u_a = t(t.u_b) + n_a(n_b.u_b)
+       !     = u_b - n_b(n_b.u_b) + n_a(n_b.u_b)
+       !     = (I + (n_a-n_b)n_b^T)u_b == B u_b
+       forall(gi=1:face_ngi(X,face))
+          Btmp(:,:,gi)=outer_product(n1(:,gi)-n2(:,gi),n2(:,gi))
+          forall(i=1:size(Btmp,1)) Btmp(i,i,gi)=Btmp(i,i,gi)+1.0
+          B(:,:,gi)=matmul(J_scaled(:,:,gi),Btmp(:,:,gi))
+       end forall
+    case (VECTOR_UPWIND_SPHERE)
+       if(mesh_dim(X).ne.2) then
+          FLExit('Assumes 2d surface. Surface of the sphere, in fact.')
+       end if
+       !Get element normal
+       up_gi = -ele_val_at_quad(down,ele)
+       call get_up_gi(X,ele,up_gi)
+       up_gi2 = -ele_val_at_quad(down,ele_2)
+       call get_up_gi(X,ele_2,up_gi2)
+       ! Form 'bending' tensor
+       ! Coordinate system is:
+       ! t1: normal to local normal and (0,0,1)
+       ! t2: normal to t1 and local normal
+       ! u2 = t21(t21.u2) + t22(t22.u2)
+       ! u1 = t11(t21.u2) + t12(t22.u2)
+       !    = (t11 t21^T + t12 t22^T)u2
+
+       pole_axis = 0.
+       pole_axis(3,:) = 1.
+       btmp = 0.
+       do gi = 1, face_ngi(X,face)
+          t11(:,gi) = cross_product(pole_axis(:,gi),up_gi(:,1))
+          t11(:,gi) = t11(:,gi)/sqrt(sum(t11(:,gi)**2))
+          t12(:,gi) = cross_product(up_gi(:,1),t11(:,gi))
+          t21(:,gi) = cross_product(pole_axis(:,gi),up_gi2(:,1))
+          t21(:,gi) = t21(:,gi)/sqrt(sum(t21(:,gi)**2))
+          t22(:,gi) = cross_product(up_gi2(:,1),t21(:,gi))
+          
+          forall(i=1:3,k=1:3)
+             Btmp(i,k,gi) = t11(i,gi)*t21(k,gi) + t12(i,gi)*t22(k,gi)
+          end forall
+          B(:,:,gi)=matmul(J_scaled(:,:,gi),Btmp(:,:,gi))
+       end do
+    case default
+       FLAbort('Unknown vector upwinding option')
+    end select
+
+    !----------------------------------------------------------------------
+    ! Construct bilinear forms.
+    !----------------------------------------------------------------------
+       
+    ! Calculate outflow boundary integral.
+    ! can anyone think of a way of optimising this more to avoid
+    ! superfluous operations (i.e. multiplying things by 0 or 1)?
+
+    ! first the integral around the inside of the element
+    ! (this is the flux *out* of the element)
+    inner_advection_integral = -income*u_nl_q_dotn_l
+    nnAdvection_out=shape_shape_tensor(U_shape, U_shape,  &
+         &inner_advection_integral*U_shape%quadrature%weight, J_scaled)
+    
+    ! now the integral around the outside of the element
+    ! (this is the flux *in* to the element)
+    outer_advection_integral = income*u_nl_q_dotn_l
+    nnAdvection_in=shape_shape_tensor(U_shape, U_shape_2, &
+         &outer_advection_integral*U_shape%quadrature%weight, B)
+       
+    !----------------------------------------------------------------------
+    ! Perform global assembly.
+    !----------------------------------------------------------------------
+
+    ! Insert advection in matrix.
+
+    ! Outflow boundary integral.
+    do dim1 = 1, dim
+       do dim2 = 1, X%dim
+          call addto(A, dim1, dim2, U_face, U_face, nnAdvection_out(dim1,dim2,:,:))
+       end do
+    end do
+
+    ! Inflow boundary integral.
+    do dim1 = 1, dim
+       do dim2 = 1, X%dim
+          call addto(A, dim1, dim2, U_face, U_face_2, nnAdvection_in(dim1,dim2,:,:))
+        end do
+    end do
+
+    contains
+
+      ! copy of outer_product in vector tools, but now as function
+      pure function outer_product(x, y) result (M)
+        !!< Give two column vectors x, y
+        !!< compute the matrix xy*
+
+        real, dimension(:), intent(in) :: x, y
+        real, dimension(size(x), size(y)) :: M
+        integer :: i, j
+
+        forall (i=1:size(x))
+           forall (j=1:size(y))
+              M(i, j) = x(i) * y(j)
+           end forall
+        end forall
+      
+      end function outer_product
+    
+  end subroutine construct_vector_adv_interface_dg
+
+end module advection_local_DG

=== added file 'assemble/Bubble_tools.F90'
--- assemble/Bubble_tools.F90	1970-01-01 00:00:00 +0000
+++ assemble/Bubble_tools.F90	2012-10-23 15:04:25 +0000
@@ -0,0 +1,485 @@
+!    Copyright (C) 2009 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 "fdebug.h"
+  
+  module bubble_tools
+    use fields
+    use field_options
+    use fields_allocates
+    use vector_tools
+    use sparse_matrices_fields
+    use solvers
+    use quadrature
+    use sparsity_patterns_meshes
+    use element_numbering
+    use state_module
+    use shape_functions
+    use quadrature
+    use shape_functions
+    use fldebug_parameters
+    use ieee_arithmetic, only: ieee_quiet_nan, ieee_value
+    use spud
+    use vtk_interfaces
+    use adjacency_lists
+    use reference_counting
+    implicit none
+
+    private
+    public get_lumped_mass_p2b,&
+         & project_to_p2b_lumped, get_p2b_lumped_mass_quadrature, &
+         & bubble_field_to_vtk
+
+#include "../femtools/Reference_count_interface_mesh_type.F90"
+
+  contains
+
+    !! This is special cased because it is only used for mass lumping
+    subroutine get_p2b_lumped_mass_quadrature(quad)
+      type(quadrature_type), intent(inout) :: quad
+      !
+      real :: wv = 1.0/20., we = 2.0/15., wg = 9./20.
+
+      call allocate(quad, 3, 7, 3)
+      quad%dim = 2
+      quad%degree = 3
+      quad%weight = 0.5*(/wv,we,wv,we,we,wv,wg/)
+      quad%l(:,1) = (/1.0,0.5,0.0,0.5,0.0,0.,1./3./)
+      quad%l(:,2) = (/0.0,0.5,1.0,0.0,0.5,0.,1./3./)
+      quad%l(:,3) = (/0.0,0.0,0.0,0.5,0.5,1.,1./3./)
+      quad%family = 3!FAMILY_SIMPSONS
+
+    end subroutine get_p2b_lumped_mass_quadrature
+
+    !! This is special cased because of the special properties of the 
+    !! p2b lumped mass (it is still 3rd order, and positive).
+    subroutine get_lumped_mass_p2b(state,mass,p2b_field)
+      type(csr_matrix), intent(inout) :: mass
+      type(state_type), intent(inout) :: state
+      type(scalar_field), intent(in) :: p2b_field
+      !
+      integer :: ele
+      type(vector_field), pointer :: X
+      real, dimension(7) :: N_vals
+      type(quadrature_type) :: Simpsons_quad
+      type(element_type) :: p2b_lumped_shape, X_lumped_shape
+
+      X=>extract_vector_field(state, "Coordinate")
+
+      if(p2b_field%mesh%shape%dim.ne.2) then
+         FLAbort('Only works for 2d meshes')
+      end if
+      if(p2b_field%mesh%shape%loc.ne.7) then
+         FLAbort('Expected p2 bubble mesh')
+      end if
+
+      call get_p2b_lumped_mass_quadrature(Simpsons_quad)
+      p2b_lumped_shape = make_element_shape_from_element(p2b_field%mesh%shape, &
+           quad=Simpsons_quad)
+      X_lumped_shape = make_element_shape_from_element(X%mesh%shape, &
+           quad=Simpsons_quad)
+
+      call zero(mass)
+
+      do ele = 1, ele_count(X)
+         call get_lumped_mass_p2b_ele(mass,p2b_field,p2b_lumped_shape,&
+              X,X_lumped_shape,&
+              ele)
+      end do
+      
+      call deallocate(p2b_lumped_shape)
+      call deallocate(X_lumped_shape)
+      call deallocate(Simpsons_quad)
+
+    end subroutine get_lumped_mass_p2b
+
+    subroutine get_lumped_mass_p2b_ele(mass,p2b_field,p2b_lumped_shape,&
+         X,X_lumped_shape,ele)
+      type(csr_matrix), intent(inout) :: mass
+      type(scalar_field), intent(in) :: p2b_field
+      type(vector_field), intent(in) :: X
+      type(element_type), intent(in) :: p2b_lumped_shape,X_lumped_shape
+      integer, intent(in) :: ele
+      !
+      real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J
+      real, dimension(ele_ngi(X,ele)) :: detwei
+      real, dimension(X_lumped_shape%ngi) :: detwei_l
+      real, dimension(ele_loc(p2b_field,ele),ele_loc(p2b_field,ele))&
+           :: l_mass_mat
+      real :: Area
+      integer :: loc
+      real, dimension(p2b_lumped_shape%ngi) :: weight_vals
+      call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), J, detwei)
+      Area = sum(detwei)
+      call compute_jacobian(ele_val(X,ele), X_lumped_shape, J, detwei_l)
+      assert(abs(Area-sum(detwei_l))<1.0e-10)
+      l_mass_mat = shape_shape(p2b_lumped_shape,p2b_lumped_shape,&
+           detwei_l)
+      call addto(mass,ele_nodes(p2b_field,ele),&
+           ele_nodes(p2b_field,ele),l_mass_mat)
+
+    end subroutine get_lumped_mass_p2b_ele
+
+    !! This is special cased because of the special properties of the 
+    !! p2b lumped mass (it is still 3rd order, and positivity preserving)
+    subroutine project_to_p2b_lumped(state,field,field_projected)
+      type(state_type), intent(inout) :: state
+      type(scalar_field), intent(in) :: field
+      type(scalar_field), intent(inout) :: field_projected
+      !
+      type(scalar_field) :: rhs
+      type(vector_field), pointer :: X
+      type(csr_sparsity), pointer :: mass_sparsity
+      type(csr_matrix) :: lumped_mass
+      integer :: ele
+      type(quadrature_type) :: Simpsons_quad
+      type(element_type) :: X_lumped_shape,&
+           & field_projected_lumped_shape, field_lumped_shape
+
+      if(field_projected%mesh%shape%dim.ne.2) then
+         FLAbort('Only works for 2d meshes')
+      end if
+      if(field_projected%mesh%shape%loc.ne.7) then
+         FLAbort('Expected p2 bubble mesh')
+      end if
+      X=>extract_vector_field(state, "Coordinate")
+
+      call get_p2b_lumped_mass_quadrature(Simpsons_quad)
+      field_projected_lumped_shape = make_element_shape_from_element( &
+           field_projected%mesh%shape,quad=Simpsons_quad)
+      field_lumped_shape = make_element_shape_from_element( &
+           field%mesh%shape,quad=Simpsons_quad)
+      X_lumped_shape = make_element_shape_from_element(X%mesh%shape, &
+           quad=Simpsons_quad)
+
+      mass_sparsity => get_csr_sparsity_firstorder(state, &
+           field_projected%mesh, field_projected%mesh)
+      call allocate(lumped_mass,mass_sparsity)
+      call zero(lumped_mass)
+      call allocate(rhs,field_projected%mesh,"ProjectionRHS")
+      call zero(rhs)
+
+      do ele = 1, ele_count(X)
+         call project_to_p2b_lumped_ele(rhs,lumped_mass,&
+              X,X_lumped_shape,field_lumped_shape,&
+              field_projected_lumped_shape,field,ele)
+      end do
+      call petsc_solve(field_projected,lumped_mass,rhs)
+      ewrite (2,*) maxval(abs(field_projected%val))
+      call deallocate(rhs)
+      call deallocate(lumped_mass)
+    end subroutine project_to_p2b_lumped
+
+    subroutine project_to_p2b_lumped_ele(rhs,lumped_mass,&
+         X,X_lumped_shape,field_lumped_shape,&
+         field_projected_lumped_shape,field,ele)
+      type(vector_field), intent(in) :: X
+      type(scalar_field), intent(in) :: field
+      type(scalar_field), intent(inout) :: rhs
+      type(csr_matrix), intent(inout) :: lumped_mass
+      integer, intent(in) :: ele
+      type(element_type), intent(in) :: X_lumped_shape,&
+           field_projected_lumped_shape, field_lumped_shape
+      !
+      real, dimension(ele_loc(rhs,ele)) :: l_rhs,field_gi      
+      real, dimension(ele_loc(rhs,ele),ele_loc(rhs,ele)) :: &
+           & l_mass
+      real, dimension(ele_loc(field,ele)) :: field_vals
+      integer :: loc
+      real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J
+      real, dimension(X_lumped_shape%ngi) :: detwei
+
+      call compute_jacobian(ele_val(X,ele), X_lumped_shape, J, detwei)
+
+      !Values of field at field DOFs
+      field_vals = ele_val(field,ele)
+      assert(size(field_gi)==size(detwei))
+      field_gi = matmul(transpose(field_lumped_shape%n),field_vals)
+      l_rhs = shape_rhs(field_projected_lumped_shape,&
+           field_gi*detwei)
+      l_mass = shape_shape(field_projected_lumped_shape,&
+           field_projected_lumped_shape,detwei)
+      call addto(rhs,ele_nodes(rhs,ele),l_rhs)
+      call addto(lumped_mass,ele_nodes(rhs,ele),&
+           ele_nodes(rhs,ele),l_mass)
+    end subroutine project_to_p2b_lumped_ele
+
+    subroutine bubble_field_to_vtk(state,sfields,filename,dump_num)
+      type(state_type), intent(inout) :: state
+      type(scalar_field), dimension(:), intent(in) :: sfields
+      character(len=*), intent(in) :: filename ! Base filename
+      integer, intent(in) :: dump_num
+      !
+      type(mesh_type), pointer :: lQuadMesh
+      type(mesh_type), target :: QuadMesh
+      type(vector_field) :: QuadX
+      type(scalar_field), dimension(size(sfields)) :: QuadSF
+      type(vector_field), pointer :: X
+      integer :: stat, i
+      type(mesh_type), pointer :: vertex_mesh
+
+      ewrite(1,*) 'subroutine bubble_field_to_vtk(state,sfield,filename,dump_num)'
+      do i = 1, size(sfields)
+         if(sfields(i)%mesh%shape%numbering%type.ne.ELEMENT_BUBBLE) then
+            FLExit('Only works for bubble functions.')
+         end if
+      end do
+      X => extract_vector_field(state,"Coordinate")
+      call find_linear_parent_mesh(state, X%mesh, vertex_mesh)
+
+      !Construct a quadrilateral mesh from the triangular one.
+      lQuadMesh => extract_mesh(state,"QuadMesh",stat)
+      if(stat/=0) then
+         call quadmesh_from_trimesh(QuadMesh,vertex_mesh,state)
+         lQuadMesh => QuadMesh
+      end if
+      !Construct a coordinate field on the quadrilateral mesh
+      call allocate(QuadX,X%dim,lQuadMesh,"QuadCoordinate")
+      call zero(QuadX)
+      call quadX_from_bubbleX(QuadX,X)
+      !Remap the bubble field to Q1 on the quad mesh
+      do i = 1, size(sfields)
+         call allocate(QuadSF(i),lQuadMesh,"Quad"//trim(sfields(i)%name))
+         call quadSF_from_bubbleSF(QuadSF(i),sfields(i))
+      end do
+      call vtk_write_fields(filename, dump_num, quadX, lQuadMesh, &
+           QuadSF)
+      call deallocate(QuadX)
+      do i = 1, size(sfields)
+         call deallocate(QuadSF(i))
+      end do
+      
+      ewrite(1,*) 'END subroutine bubble_field_to_vtk(state,sfield,filename,dump_num)'
+      contains 
+
+        subroutine quadmesh_from_trimesh(QuadMesh,Xmesh,state)
+          type(mesh_type), intent(inout) :: QuadMesh
+          type(mesh_type), intent(in) :: Xmesh
+          type(state_type), intent(inout) :: state
+          !
+          type(csr_sparsity) :: NNList, EEList
+          type(csr_matrix) :: edge_matrix
+          type(quadrature_type) :: quad
+          integer :: tri_edges, tri_nodes, tri_elements
+          integer :: ele, ni, nodecount, node, qele, qele_loc,ele2
+          integer, pointer, dimension(:) :: neigh, X_ele
+
+          tri_elements = Xmesh%elements
+          quadmesh%elements = 3*tri_elements
+          allocate(quadmesh%ndglno(quadmesh%elements*4))
+          quadmesh%ndglno=-666
+          call MakeLists_Mesh(Xmesh,NNList=NNlist,EEList=EEList)
+          call allocate(edge_matrix,EEList,type=CSR_INTEGER)
+          call zero(edge_matrix)
+          tri_edges=size(NNList%colm)/2
+          tri_nodes=Xmesh%nodes
+          quadmesh%nodes = 3*tri_elements+tri_edges+Xmesh%nodes
+
+          nodecount = Xmesh%nodes
+          do ele = 1, ele_count(Xmesh)             
+             !! Element numbering in the three quads
+             !! ele_q = 3*(ele-1)+vnode where vnode is the local
+             !! node number of the node at the vertex
+             !!Local numbering in the three quads
+             !!triangle vertex node first, then edge node
+             !!on edge joining to vertex with next number in modulo arithmetic
+             !!then triangle centre, then other edge
+
+             !First do the vertex nodes (same numbering as Xmesh)
+             X_ele => ele_nodes(Xmesh,ele)
+             assert(size(X_ele)==3)
+             do qele_loc = 1,3
+                qele = 3*(ele-1)+qele_loc
+                quadmesh%ndglno(4*(qele-1)+1) = X_ele(qele_loc)
+             end do
+             
+             !Then the edge nodes (only do if ele<ele2 or ele2<0)
+             neigh => ele_neigh(Xmesh,ele)
+             assert(size(neigh)==3)
+             do ni = 1, size(neigh)
+                ele2 = neigh(ni)
+                if(ele2.le.0 .or. ele<ele2) then
+                   nodecount = nodecount + 1
+                   node = nodecount
+                   if(ele2>0) then
+                      call set(edge_matrix,ele,ele2,node)
+                      call set(edge_matrix,ele2,ele,node)
+                   end if
+                else
+                   node = val(edge_matrix,ele,ele2)
+                end if
+                assert(node>0)
+                select case(ni)
+                case (1) 
+                   qele = 3*(ele-1)+2
+                   quadmesh%ndglno(4*(qele-1)+2)=node
+                   qele = 3*(ele-1)+3
+                   quadmesh%ndglno(4*(qele-1)+3)=node
+                case (2)
+                   qele = 3*(ele-1)+3
+                   quadmesh%ndglno(4*(qele-1)+2)=node
+                   qele = 3*(ele-1)+1
+                   quadmesh%ndglno(4*(qele-1)+3)=node
+                case (3)
+                   qele = 3*(ele-1)+1
+                   quadmesh%ndglno(4*(qele-1)+2)=node
+                   qele = 3*(ele-1)+2
+                   quadmesh%ndglno(4*(qele-1)+3)=node
+                end select
+             end do
+
+             !Then the triangle centre nodes
+             do qele = 3*(ele-1)+1,3*(ele-1)+3
+                nodecount = nodecount + 1
+                quadmesh%ndglno(4*(qele-1)+4)=nodecount
+             end do
+          end do
+          assert(all(quadmesh%ndglno>0))
+          assert(nodecount==quadmesh%nodes)
+
+          quadmesh%wrapped=.false.
+          quad = make_quadrature(4,2,degree=Xmesh%shape%quadrature%degree,&
+               family=Xmesh%shape%quadrature%family)
+          quadmesh%shape = make_element_shape(4,2,1,quad)
+          quadmesh%name = "QuadMesh"
+          quadmesh%continuity=0
+          quadmesh%periodic=.false.
+
+          allocate(quadmesh%adj_lists)
+          nullify(quadmesh%region_ids)
+          nullify(quadmesh%subdomain_mesh)
+          nullify(quadmesh%refcount) ! Hack for gfortran component initialisation
+          !                         bug.
+
+          call addref(quadmesh)
+          
+          !need to insert into state and deallocate
+          call insert(state,quadmesh,"QuadMesh")
+          call deallocate(quadmesh)
+          call deallocate(EEList)
+          call deallocate(NNList)
+        end subroutine quadmesh_from_trimesh
+
+        subroutine quadX_from_bubbleX(QuadX,X)
+          type(vector_field), intent(inout) :: QuadX
+          type(vector_field), intent(in) :: X
+          !
+          real, dimension(X%dim,ele_loc(X,1)) :: X_vals
+          real, dimension(X%dim,ele_loc(QuadX,1)) :: QX_vals
+          integer :: qele,dim1,ele
+          
+          real, dimension(X%dim,3) :: X_lin
+
+          do ele = 1, ele_count(X)
+             X_vals = ele_val(X,ele)
+
+             !!Need to evaluate X at bubble node locations
+             select case(ele_loc(X,1))
+                case(3)
+                   X_lin = X_vals
+                case(6)
+                   X_lin(:,1) = X_vals(:,1)
+                   X_lin(:,2) = X_vals(:,3)
+                   X_lin(:,3) = X_vals(:,6)
+                case(10)
+                   X_lin(:,1) = X_vals(:,1)
+                   X_lin(:,2) = X_vals(:,4)
+                   X_lin(:,3) = X_vals(:,10)
+                case default
+                   FLExit('Dont know how to handle coordinate mesh')
+             end select
+
+             !! Quadrilateral #1
+             do dim1 = 1, X%dim
+                QX_vals(dim1,:) = &
+                     &(/X_lin(dim1,1),0.5*(X_lin(dim1,1)+X_lin(dim1,2)),&
+                     &0.5*(X_lin(dim1,1)+X_lin(dim1,3)),&
+                     &sum(X_lin(dim1,:))/3/)
+             end do
+             qele = (ele-1)*3+1
+             call set(QuadX,ele_nodes(QuadX,qele),QX_vals)
+             !! Quadrilateral #2
+             do dim1 = 1, X%dim
+                QX_vals(dim1,:) = &
+                     &(/X_lin(dim1,2),0.5*(X_lin(dim1,2)+X_lin(dim1,3)),&
+                     &0.5*(X_lin(dim1,1)+X_lin(dim1,2)),&
+                     &sum(X_lin(dim1,:))/3/)
+             end do
+             qele = (ele-1)*3+2
+             call set(QuadX,ele_nodes(QuadX,qele),QX_vals)                
+             !! Quadrilateral #3
+             do dim1 = 1, X%dim
+                QX_vals(dim1,:) = &
+                     &(/X_lin(dim1,3),0.5*(X_lin(dim1,3)+X_lin(dim1,1)),&
+                     &0.5*(X_lin(dim1,2)+X_lin(dim1,3)),&
+                     &sum(X_lin(dim1,:))/3/)
+             end do
+             qele = (ele-1)*3+3
+             call set(QuadX,ele_nodes(QuadX,qele),QX_vals)                
+          end do
+
+        end subroutine quadX_from_bubbleX
+
+        subroutine quadSF_from_bubbleSF(QuadSF,sfield)
+          type(scalar_field), intent(inout) :: QuadSF
+          type(scalar_field), intent(in) :: sfield
+          
+          !
+          real, dimension(ele_loc(sfield,1)) :: s_vals
+          real, dimension(ele_loc(QuadSF,1)) :: qs_vals
+          integer :: qele,ele
+          real, dimension(7) :: N_vals
+          
+          !Basis functions evaluated at bubble node.
+          N_vals = eval_shape(ele_shape(sfield,1), (/1.0/3.0,1.0/3.0,1.0/3.0/))
+
+          do ele = 1, ele_count(sfield)
+             s_vals = ele_val(sfield,ele)
+             assert(size(s_vals)==7)
+             s_vals(7) = s_vals(7)*N_vals(7)
+             s_vals(7) = s_vals(7) + sum(s_vals(1:6)*N_vals(1:6))
+             !! Quadrilateral #1
+             QS_vals(:) = &
+                  &(/s_vals(1),s_vals(2),s_vals(4),s_vals(7)/)
+             qele = (ele-1)*3+1
+             call set(QuadSF,ele_nodes(QuadSF,qele),QS_vals)
+             !! Quadrilateral #2
+             QS_vals(:) = &
+                  &(/s_vals(3),s_vals(5),s_vals(2),s_vals(7)/)
+             qele = (ele-1)*3+2
+             call set(QuadSF,ele_nodes(QuadSF,qele),QS_vals)
+             !! Quadrilateral #3
+             QS_vals(:) = &
+                  &(/s_vals(6),s_vals(4),s_vals(5),s_vals(7)/)
+             qele = (ele-1)*3+3
+             call set(QuadSF,ele_nodes(QuadSF,qele),QS_vals)
+          end do
+        end subroutine quadSF_from_bubbleSF
+      end subroutine bubble_field_to_vtk
+
+#include "../femtools/Reference_count_mesh_type.F90"
+      
+  end module bubble_tools

=== modified file 'assemble/Diagnostic_fields_wrapper.F90'
--- assemble/Diagnostic_fields_wrapper.F90	2012-09-18 08:03:27 +0000
+++ assemble/Diagnostic_fields_wrapper.F90	2012-10-23 15:04:25 +0000
@@ -79,7 +79,7 @@
     ! An array of submaterials of the current phase in state(istate).
     type(state_type), dimension(:), pointer :: submaterials
     
-    ewrite(1, *) "In calculate_diagnostic_variables"
+    ewrite(1, *) "In calculate_diagnostic_variables (Diagnostic_fields_wrapper)"
  
     do i = 1, size(state)
 

=== modified file 'assemble/Hybridized_Helmholtz.F90'
--- assemble/Hybridized_Helmholtz.F90	2011-07-22 12:39:08 +0000
+++ assemble/Hybridized_Helmholtz.F90	2012-10-23 15:04:25 +0000
@@ -24,7 +24,6 @@
   !    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"
 
@@ -44,18 +43,213 @@
     use diagnostic_fields_wrapper
     use assemble_cmc
     use FUtils, only : real_vector, real_matrix
-    use global_parameters, only: option_path_len
+    use global_parameters, only: option_path_len, PYTHON_FUNC_LEN
     use vector_tools, only: solve
-    use manifold_projections
-    implicit none
-
-contains
+    use manifold_tools
+    implicit none
+
+    !Variables belonging to the Newton solver
+    logical, private :: newton_initialised = .false.
+    !Cached Newton solver matrices
+    type(csr_matrix), private :: Newton_lambda_mat
+    type(block_csr_matrix), private :: Newton_continuity_mat
+    type(real_matrix), pointer, dimension(:), private &
+         &:: Newton_local_solver_cache=>null()
+    type(real_matrix), pointer, dimension(:), private &
+         &:: Newton_local_solver_rhs_cache=>null()
+
+contains 
+
+  subroutine solve_hybridised_timestep_residual(state,newU,newD,&
+       UResidual, DResidual)
+
+    ! Subroutine to apply one Newton iteration to hybridised shallow
+    ! water equations
+    ! Written to cache all required matrices
+    implicit none
+    type(state_type), intent(inout) :: state
+    type(vector_field), intent(inout) :: newU !U at next timestep
+    type(scalar_field), intent(inout) :: newD !D at next timestep
+    type(vector_field), intent(in), optional :: UResidual
+    type(scalar_field), intent(in), optional :: DResidual
+    !
+    type(vector_field), pointer :: X, U, down, U_cart, U_old
+    type(scalar_field), pointer :: D,f, D_old
+    type(scalar_field) :: lambda, D_res, D_res2
+    type(vector_field) :: U_res, U_res2
+    type(scalar_field), target :: lambda_rhs, u_cpt
+    type(csr_sparsity) :: lambda_sparsity, continuity_sparsity
+    type(csr_matrix) :: continuity_block_mat
+    type(mesh_type), pointer :: lambda_mesh
+    real :: D0, dt, g, theta, tolerance
+    integer :: ele,i1, stat, dim1, dloc, uloc,mdim,n_constraints
+    logical :: have_constraint
+    character(len=OPTION_PATH_LEN) :: constraint_option_string
+
+    ewrite(1,*) 'solve_hybridised_timestep_residual(state)'
+
+    ewrite(1,*) 'TIMESTEPPING LAMBDA'
+
+    !get parameters
+    call get_option("/physical_parameters/gravity/magnitude", g)
+    call get_option("/timestepping/theta",theta)
+    call get_option("/material_phase::Fluid/scalar_field::LayerThickness/&
+         &prognostic/mean_layer_thickness",D0)
+    call get_option("/timestepping/timestep", dt)
+
+    !Pull the fields out of state
+    D=>extract_scalar_field(state, "LayerThickness")
+    f=>extract_scalar_field(state, "Coriolis")
+    U=>extract_vector_field(state, "LocalVelocity")
+    X=>extract_vector_field(state, "Coordinate")
+    down=>extract_vector_field(state, "GravityDirection")
+    U_cart => extract_vector_field(state, "Velocity")
+    U_old => extract_vector_field(state, "OldLocalVelocity")
+    D_old => extract_scalar_field(state, "OldLayerThickness")
+
+    !Allocate local field variables
+    lambda_mesh=>extract_mesh(state, "VelocityMeshTrace")
+    call allocate(lambda,lambda_mesh,name="LagrangeMultiplier")
+    call allocate(lambda_rhs,lambda%mesh,"LambdaRHS")
+    call zero(lambda_rhs)
+    call allocate(U_res,U%dim,U%mesh,"VelocityResidual")
+    call allocate(D_res,D%mesh,"LayerThicknessResidual")
+    call allocate(U_res2,U%dim,U%mesh,"VelocityResidual2")
+    call allocate(D_res2,D%mesh,"LayerThicknessResidual2")
+
+    if(.not.newton_initialised) then
+       !construct/extract sparsities
+       lambda_sparsity=get_csr_sparsity_firstorder(&
+            &state, lambda%mesh, lambda&
+            &%mesh)
+       continuity_sparsity=get_csr_sparsity_firstorder(state, u%mesh,&
+            & lambda%mesh)
+
+       !allocate matrices
+       call allocate(Newton_lambda_mat,lambda_sparsity)
+       call zero(Newton_lambda_mat)
+       call allocate(Newton_continuity_mat,continuity_sparsity,(/U%dim,1/))
+       call zero(Newton_continuity_mat)
+
+       mdim = mesh_dim(U)
+       have_constraint = &
+            &have_option(trim(U%mesh%option_path)//"/from_mesh/constraint_type")
+
+       allocate(newton_local_solver_cache(ele_count(U)))
+       allocate(newton_local_solver_rhs_cache(ele_count(U)))
+       do ele = 1, ele_count(D)
+          uloc = ele_loc(U,ele)
+          dloc = ele_loc(d,ele)
+          n_constraints = 0
+          if(have_constraint) then
+             n_constraints = ele_n_constraints(U,ele)
+             allocate(newton_local_solver_cache(ele)%ptr(&
+                  mdim*uloc+dloc+n_constraints,&
+                  mdim*uloc+dloc+n_constraints))
+             allocate(newton_local_solver_rhs_cache(ele)%ptr(&
+                  mdim*uloc+dloc,mdim*uloc+dloc))
+          end if
+       end do
+
+       !Assemble matrices
+       do ele = 1, ele_count(D)
+          call assemble_newton_solver_ele(D,f,U,X,down,lambda_rhs,ele, &
+               &g,dt,theta,D0,&
+               &lambda_mat=Newton_lambda_mat,&
+               &continuity_mat=Newton_continuity_mat,&
+               &Local_solver_matrix=Newton_local_solver_cache(ele)%ptr,&
+               &Local_solver_rhs=Newton_local_solver_rhs_cache(ele)%ptr)
+       end do
+       newton_initialised = .true.
+    end if
+
+    if(present(DResidual).and.present(UResidual).and..not.&
+         have_option('/material_phase::Fluid/vector_field::Velocity/prognostic/wave_equation/fully_coupled/linear_debug')) then
+       !Set residuals from nonlinear input
+       call set(U_res,UResidual)
+       call set(D_res,DResidual)
+    else
+       !Compute residuals for linear equation
+       do ele = 1, ele_count(D)
+          call get_linear_residuals_ele(U_res,D_res,&
+               U_old,D_old,newU,newD,&
+               newton_local_solver_cache(ele)%ptr,&
+               newton_local_solver_rhs_cache(ele)%ptr,ele)
+       end do
+       ewrite(2,*) 'cjc U_res calc', sum(newU%val), maxval(abs(U_res%val))
+       ewrite(2,*) 'cjc D_res calc', sum(newD%val), maxval(abs(D_res%val))
+    end if
+
+    do ele = 1, ele_count(D)
+       call local_solve_residuals_ele(U_res,D_res,&
+            newton_local_solver_cache(ele)%ptr,ele)
+    end do
+
+    call mult_t(lambda_rhs,Newton_continuity_mat,U_res)
+    call scale(lambda_rhs,-1.0)
+    ewrite(2,*) 'LAMBDARHS', maxval(abs(lambda_rhs%val))
+
+    call zero(lambda)
+    !Solve the equations
+    call petsc_solve(lambda,newton_lambda_mat,lambda_rhs,&
+         option_path=trim(U_cart%mesh%option_path)//&
+         &"/from_mesh/constraint_type")
+    ewrite(2,*) 'LAMBDA', maxval(abs(lambda%val))
+
+    !Update new U and new D from lambda
+    !Compute residuals 
+    if(present(DResidual).and.present(UResidual).and..not.&
+         have_option('/material_phase::Fluid/vector_field::Velocity/prognostic/wave_equation/fully_coupled/linear_debug')) then
+       call set(U_res,UResidual)
+       call set(D_res,DResidual)
+    else
+       do ele = 1, ele_count(D)
+          call get_linear_residuals_ele(U_res,D_res,&
+               U_old,D_old,newU,newD,&
+               newton_local_solver_cache(ele)%ptr,&
+               newton_local_solver_rhs_cache(ele)%ptr,ele)
+       end do
+    end if
+    ewrite(2,*) 'cjc U_res recalc', sum(newU%val), maxval(abs(U_res%val))
+    ewrite(2,*) 'cjc D_res recalc', sum(newD%val), maxval(abs(D_res%val))
+       
+    call mult(U_res2,Newton_continuity_mat,lambda)
+    call addto(U_res,U_res2)
+    do ele = 1, ele_count(U)
+       call local_solve_residuals_ele(U_res,D_res,&
+            newton_local_solver_cache(ele)%ptr,ele)
+    end do
+
+    !Negative scaling occurs here.
+    call scale(U_res,-1.0)
+    call scale(D_res,-1.0)
+
+    ewrite(2,*) 'Delta U', maxval(abs(U_res%val))
+    ewrite(2,*) 'Delta D', maxval(abs(D_res%val))
+    
+    call addto(newU,U_res)
+    call addto(newD,D_res)
+
+    !Deallocate local variables
+    call deallocate(lambda_rhs)
+    call deallocate(lambda)
+    call deallocate(U_res)
+    call deallocate(D_res)
+    call deallocate(U_res2)
+    call deallocate(D_res2)
+
+    ewrite(1,*) 'END solve_hybridised_timestep_residual(state)'
+  
+  end subroutine solve_hybridised_timestep_residual
+
   subroutine solve_hybridized_helmholtz(state,D_rhs,U_Rhs,&
        &D_out,U_out,&
+       &dt_in,theta_in,&
        &compute_cartesian,&
-       &check_continuity,output_dense,&
-       &projection,poisson,u_rhs_local,&
-       &solver_option_path)
+       &output_dense,&
+       &projection,poisson,&
+       &u_rhs_local)
+
     ! Subroutine to solve hybridized helmholtz equation
     ! If D_rhs (scalar pressure field) is present, then solve:
     ! <w,u> + <w,fu^\perp> - g <div w,d> + <<[w],d>> = <w,U_rhs>
@@ -63,7 +257,7 @@
     ! <<\gamma, [u]>> = 0
     ! (i.e. for unit testing)
     ! otherwise solve:
-    ! <w,u> + dt*theta*<w,fu^\perp> - dt*theta*g <div w,d> + <<[w],d>> = 
+    ! <w,u> + dt*theta*<w,fu^\perp> - dt*theta*g <div w,d> + <<[w],l>> = 
     ! -dt*<w f(u^n)^\perp> + dt*g*<div w, d^n>
     ! <\phi,\eta> +  dt*theta*<\phi,div u> = <\ph
     ! <<\gamma, [u]>> = 0
@@ -76,10 +270,10 @@
     type(vector_field), intent(inout), optional :: U_rhs
     type(scalar_field), intent(inout), optional :: D_out
     type(vector_field), intent(inout), optional :: U_out
+    real, intent(in), optional :: theta_in,dt_in
     logical, intent(in), optional :: compute_cartesian, &
-         &check_continuity,output_dense, projection,poisson
+         &output_dense, projection,poisson
     logical, intent(in), optional :: u_rhs_local !means u_rhs is in local coords
-    character(len=OPTION_PATH_LEN), intent(in), optional :: solver_option_path
     !
     type(vector_field), pointer :: X, U, down, U_cart
     type(scalar_field), pointer :: D,f, lambda_nc
@@ -89,21 +283,36 @@
     type(csr_matrix) :: lambda_mat, continuity_block_mat,continuity_block_mat1
     type(block_csr_matrix) :: continuity_mat
     type(mesh_type), pointer :: lambda_mesh
-    real :: D0, dt, g, theta
+    real :: D0, dt, g, theta, tolerance
     integer :: ele,i1, stat, dim1
-    logical :: l_compute_cartesian,l_check_continuity, l_output_dense
+    logical :: l_compute_cartesian, l_output_dense
     real, dimension(:,:), allocatable :: lambda_mat_dense
     character(len=OPTION_PATH_LEN) :: constraint_option_string
+    real :: u_max
+    real, dimension(:), allocatable :: weights
+    logical :: l_projection,l_poisson, pullback
 
     ewrite(1,*) '  subroutine solve_hybridized_helmholtz('
 
     l_compute_cartesian = .false.
     if(present(compute_cartesian)) l_compute_cartesian = compute_cartesian
-    if(present(check_continuity)) l_check_continuity = check_continuity
-    if(l_check_continuity) l_compute_cartesian = .true.
     l_output_dense = .false.
     if(present(output_dense)) l_output_dense = output_dense
 
+    if(present(projection)) then
+       l_projection = projection
+    else
+       l_projection = .false.
+    end if
+    if(present(poisson)) then
+       l_poisson = poisson
+    end if
+    ewrite(2,*) 'Projection = ', l_projection
+    ewrite(2,*) 'Poisson = ', l_poisson
+    if(l_poisson.and.l_projection) then
+       FLAbort('Can''t do projection and poisson')
+    end if
+
     !Pull the fields out of state
     D=>extract_scalar_field(state, "LayerThickness")
     f=>extract_scalar_field(state, "Coriolis")
@@ -134,68 +343,86 @@
     !get parameters
     call get_option("/physical_parameters/gravity/magnitude", g)
     !theta
-    call get_option("/material_phase::Fluid/scalar_field::LayerThickness/&
-         &prognostic/temporal_discretisation/theta",theta)
+
+    if(present(theta_in)) then
+       theta = theta_in
+    else
+       call get_option("/timestepping/theta",theta)
+    end if
     !D0
     call get_option("/material_phase::Fluid/scalar_field::LayerThickness/&
-         &p&
-         &rognostic/mean_layer_thickness",D0)
-    call get_option("/timestepping/timestep", dt)
-    
+         &prognostic/mean_layer_thickness",D0)
+    if(present(dt_in)) then
+       dt = dt_in
+    else
+       call get_option("/timestepping/timestep", dt)
+    end if
+
+    !compute rescaling weights for local coordinates
+    allocate(weights(element_count(X)))
+    weights = 1.0
+    !call get_weights(X,weights)
+
+    pullback = have_option('/material_phase::Fluid/scalar_field::LayerThickn&
+         &ess/prognostic/spatial_discretisation/discontinuous_galerkin/wave_&
+         &equation/pullback')
+
     !Assemble matrices
     do ele = 1, ele_count(D)
        call assemble_hybridized_helmholtz_ele(D,f,U,X,down,ele, &
-            &g,dt,theta,D0,lambda_mat=lambda_mat,&
+            &g,dt,theta,D0,pullback,weights(ele),lambda_mat=lambda_mat,&
             &lambda_rhs=lambda_rhs,D_rhs=D_rhs,U_rhs=U_rhs,&
             &continuity_mat=continuity_mat,&
-            &projection=projection,poisson=poisson,&
+            &projection=l_projection,poisson=l_poisson,&
             &u_rhs_local=u_rhs_local)
     end do
 
-    ewrite(1,*) 'LAMBDARHS', maxval(abs(lambda_rhs%val))
-
+    ewrite(2,*) 'LAMBDARHS', maxval(abs(lambda_rhs%val))
+    call zero(lambda)
     !Solve the equations
-    if(present(solver_option_path)) then
-       call petsc_solve(lambda,lambda_mat,lambda_rhs,&
-            option_path=solver_option_path)
-    else
-       call petsc_solve(lambda,lambda_mat,lambda_rhs,&
-            option_path=trim(U_cart%mesh%option_path)//&
-            &"/from_mesh/constraint_type")
-    end if
-    ewrite(1,*) 'LAMBDA', maxval(abs(lambda%val))
+    call petsc_solve(lambda,lambda_mat,lambda_rhs,&
+         option_path=trim(U_cart%mesh%option_path)//&
+         &"/from_mesh/constraint_type")
 
     !Reconstruct U and D from lambda
     do ele = 1, ele_count(D)
        call reconstruct_u_d_ele(D,f,U,X,down,ele, &
-            &g,dt,theta,D0,D_rhs=D_rhs,U_rhs=U_rhs,lambda=lambda,&
+            &g,dt,theta,D0,pullback,weights(ele),&
+            &D_rhs=D_rhs,U_rhs=U_rhs,lambda=lambda,&
             &D_out=D_out,U_out=U_out,&
-            &projection=projection,poisson=poisson,&
+            &projection=l_projection,poisson=l_poisson,&
             &u_rhs_local=u_rhs_local)
     end do
 
+    if(l_poisson.and.present(D_out)) then
+       call check_zero_level(D_out)
+    end if
+    ewrite(2,*) 'LAMBDA', maxval(abs(lambda%val))
+
     if(l_output_dense) then
        allocate(lambda_mat_dense(node_count(lambda),node_count(lambda)))
        lambda_mat_dense = dense(lambda_mat)
-       ewrite(1,*) '-----------'
+       ewrite(2,*) '-----------'
        do i1 = 1, node_count(lambda)
-          ewrite(1,*) lambda_mat_dense(i1,:)
+          ewrite(2,*) lambda_mat_dense(i1,:)
        end do
-       ewrite(1,*) '-----------'
+       ewrite(2,*) '-----------'
     end if
 
     if(l_compute_cartesian) then
        U_cart => extract_vector_field(state, "Velocity")
        if(present(U_out)) then
-          call project_local_to_cartesian(X,U_out,U_cart)
+          call project_local_to_cartesian(X,U_out,U_cart,weights=weights)
        else
-          call project_local_to_cartesian(X,U,U_cart)
+          call project_local_to_cartesian(X,U,U_cart,weights=weights)
        end if
     end if
-    if(l_check_continuity) then
-       ewrite(1,*) 'Checking continuity'
+
+    if(have_option('/geometry/mesh::VelocityMesh/check_continuity_matrix')) then
+       ewrite(2,*) 'Checking continuity'
 
        call zero(lambda_rhs)
+       u_max = 0.
        do dim1 = 1,U%dim
           if(present(U_out)) then
              u_cpt = extract_scalar_field(U_out,dim1)
@@ -204,29 +431,264 @@
           end if
           continuity_block_mat = block(continuity_mat,dim1,1)
           call mult_T_addto(lambda_rhs,continuity_block_mat,u_cpt)
-          ewrite(1,*) 'U, lambda',&
+          ewrite(2,*) 'U, lambda',&
                &maxval(u_cpt%val), maxval(abs(lambda_rhs%val))
+          u_max = u_max + maxval(abs(u_cpt%val))
        end do
-       ewrite(1,*)'JUMPS MIN:MAX',minval(lambda_rhs%val),&
-            &maxval(lambda_rhs%val)
-       assert(maxval(abs(lambda_rhs%val))<1.0e-10)
-       
-       ewrite(1,*) 'D MAXABS', maxval(abs(D%val))
+       ewrite(2,*)'JUMPS MIN:MAX',minval(lambda_rhs%val),&
+            &maxval(lambda_rhs%val), u_max
+
+       call get_option('/geometry/mesh::VelocityMesh/check_continuity_matrix/tolerance',tolerance)
+       if(maxval(abs(lambda_rhs%val))/max(1.0,u_max/3.0)>tolerance) then
+          ewrite(-1,*) 'value =', maxval(abs(lambda_rhs%val))/max(1.0,u_max/3.0)
+          FLExit('Continuity matrix tolerance failure')
+       end if
+    end if
+
+    if(have_option('/geometry/mesh::VelocityMesh/check_continuity')) then       
+       U_cart => extract_vector_field(state, "Velocity")
+       if(present(U_out)) then
+          call project_local_to_cartesian(X,U_out,U_cart,weights=weights)
+       else
+          call project_local_to_cartesian(X,U,U_cart,weights=weights)
+       end if
+       call get_option('/geometry/mesh::VelocityMesh/check_continuity/tolerance', tolerance)
        do ele = 1, ele_count(U)
-          call check_continuity_ele(U_cart,X,ele)
+          call check_continuity_ele(U_cart,X,ele,tolerance)
        end do
     end if
     
     call deallocate(lambda_mat)
+    call deallocate(continuity_mat)
     call deallocate(lambda_rhs)
     call deallocate(lambda)
+    deallocate(weights)
 
     ewrite(1,*) 'END subroutine solve_hybridized_helmholtz'
 
   end subroutine solve_hybridized_helmholtz
+
+  subroutine get_linear_residuals_ele(U_res,D_res,U,D,&
+       newU,newD,local_solver,local_solver_rhs,ele)
+    !Subroutine
+    type(vector_field), intent(inout) :: U_res,U,newU
+    type(scalar_field), intent(inout) :: D_res,D,newD
+    real, dimension(:,:), intent(in) :: local_solver, local_solver_rhs
+    integer, intent(in) :: ele
+    !
+    integer :: uloc, dloc, dim1, d_start, d_end, mdim
+    integer, dimension(mesh_dim(U)) :: U_start, U_end
+    real, dimension(mesh_dim(U)*ele_loc(U,ele)+ele_loc(D,ele))&
+         :: rhs_loc1,rhs_loc2
+    real, dimension(mesh_dim(U),ele_loc(U,ele)) :: U_val
+
+    !Calculate indices in a vector containing all the U and D dofs in
+    !element ele, First the u1 components, then the u2 components, then the
+    !D components are stored.
+    uloc = ele_loc(U_res,ele)
+    dloc = ele_loc(D_res,ele)
+    d_end = mesh_dim(U_res)*uloc+dloc
+    mdim = mesh_dim(U)
+    do dim1 = 1, mdim
+       u_start(dim1) = uloc*(dim1-1)+1
+       u_end(dim1) = uloc*dim1
+    end do
+    d_start = uloc*mdim + 1
+    d_end   = uloc*mdim+dloc
+
+    U_val = ele_val(U,ele)
+    do dim1 = 1, mesh_dim(U)
+       rhs_loc1(u_start(dim1):u_end(dim1)) = u_val(dim1,:)
+    end do
+    rhs_loc1(d_start:d_end) = ele_val(D,ele)
+    rhs_loc1 = matmul(local_solver_rhs,rhs_loc1)
+    U_val = ele_val(newU,ele)
+    do dim1 = 1, mesh_dim(U)
+       rhs_loc2(u_start(dim1):u_end(dim1)) = u_val(dim1,:)
+    end do
+    rhs_loc2(d_start:d_end) = ele_val(newD,ele)
+    rhs_loc2 = matmul(local_solver(1:d_end,1:d_end),rhs_loc2)
+    rhs_loc2 = rhs_loc2-rhs_loc1
+    do dim1 = 1, mesh_dim(U)
+       call set(U_res,dim1,ele_nodes(U,ele),&
+            rhs_loc2(u_start(dim1):u_end(dim1)))
+    end do
+    call set(D_res,ele_nodes(D,ele),&
+         rhs_loc2(d_start:d_end))
+  end subroutine get_linear_residuals_ele
+
+  subroutine local_solve_residuals_ele(U_res,D_res,&
+       local_solver_matrix,ele)
+    type(vector_field), intent(inout) :: U_res
+    type(scalar_field), intent(inout) :: D_res
+    real, dimension(:,:), intent(in) :: local_solver_matrix
+    integer, intent(in) :: ele
+    !
+    real, dimension(mesh_dim(U_res)*ele_loc(U_res,ele)+ele_loc(D_res,ele)+&
+         &ele_n_constraints(U_res,ele)) :: rhs_loc
+    real, dimension(mesh_dim(U_res),ele_loc(U_res,ele)) :: U_val
+    integer :: uloc,dloc, d_start, d_end, dim1, mdim
+    integer, dimension(mesh_dim(U_res)) :: U_start, U_end
+
+    rhs_loc = 0.
+
+    !Calculate indices in a vector containing all the U and D dofs in
+    !element ele, First the u1 components, then the u2 components, then the
+    !D components are stored.
+    uloc = ele_loc(U_res,ele)
+    dloc = ele_loc(D_res,ele)
+    d_end = mesh_dim(U_res)*uloc+dloc
+    mdim = mesh_dim(U_res)
+    do dim1 = 1, mdim
+       u_start(dim1) = uloc*(dim1-1)+1
+       u_end(dim1) = uloc*dim1
+    end do
+    d_start = uloc*mdim + 1
+    d_end   = uloc*mdim+dloc
+
+    U_val = ele_val(U_res,ele)
+    do dim1 = 1, mesh_dim(U_res)
+       rhs_loc(u_start(dim1):u_end(dim1)) = u_val(dim1,:)
+    end do
+    rhs_loc(d_start:d_end) = &
+         & ele_val(D_res,ele)
+
+    call solve(local_solver_matrix,Rhs_loc)
+    do dim1 = 1, mesh_dim(U_res)
+       call set(U_res,dim1,ele_nodes(U_res,ele),&
+            rhs_loc(u_start(dim1):u_end(dim1)))
+    end do
+    call set(D_res,ele_nodes(D_res,ele),&
+         rhs_loc(d_start:d_end))
+  end subroutine local_solve_residuals_ele
+
+  subroutine assemble_newton_solver_ele(D,f,U,X,down,lambda_rhs,ele, &
+       &g,dt,theta,D0,&
+       &lambda_mat,continuity_mat,&
+       &Local_solver_matrix,Local_solver_rhs)
+    implicit none
+    type(scalar_field), intent(in) :: D,f
+    type(scalar_field), intent(in) :: lambda_rhs
+    type(vector_field), intent(in) :: U,X,down
+    integer, intent(in) :: ele
+    real, intent(in) :: g,dt,theta,D0
+    type(csr_matrix), intent(inout) :: lambda_mat
+    type(block_csr_matrix), intent(inout) :: continuity_mat
+    real, dimension(:,:), intent(inout) ::&
+         & local_solver_matrix
+    real, dimension(:,:), intent(inout) ::&
+         & local_solver_rhs
+    !
+    real, allocatable, dimension(:,:),target :: &
+         &l_continuity_mat, l_continuity_mat2
+    real, allocatable, dimension(:,:) :: helmholtz_loc_mat
+    real, allocatable, dimension(:,:,:) :: continuity_face_mat
+    real, allocatable, dimension(:,:) :: scalar_continuity_face_mat
+    integer :: ni, face
+    integer, dimension(:), pointer :: neigh
+    type(element_type) :: U_shape
+    integer :: stat, d_start, d_end, dim1, mdim, uloc,dloc, lloc, iloc
+    integer, dimension(mesh_dim(U)) :: U_start, U_end
+    type(real_matrix), dimension(mesh_dim(U)) :: &
+         & continuity_mat_u_ptr
+    logical :: have_constraint
+    integer :: constraint_choice, n_constraints, constraints_start
+
+    !Get some sizes
+    lloc = ele_loc(lambda_rhs,ele)
+    mdim = mesh_dim(U)
+    uloc = ele_loc(U,ele)
+    dloc = ele_loc(d,ele)
+    U_shape = ele_shape(U,ele)
+
+    have_constraint = &
+         &have_option(trim(U%mesh%option_path)//"/from_mesh/constraint_type")
+    n_constraints = 0
+    if(have_constraint) then
+       n_constraints = ele_n_constraints(U,ele)
+    end if
+
+    !Calculate indices in a vector containing all the U and D dofs in
+    !element ele, First the u1 components, then the u2 components, then the
+    !D components are stored.
+    do dim1 = 1, mdim
+       u_start(dim1) = uloc*(dim1-1)+1
+       u_end(dim1) = uloc*dim1
+    end do
+    d_start = uloc*mdim + 1
+    d_end   = uloc*mdim+dloc
+
+    !Get pointers to different parts of l_continuity_mat
+    allocate(l_continuity_mat(2*uloc+dloc+n_constraints,lloc))
+    do dim1= 1,mdim
+       continuity_mat_u_ptr(dim1)%ptr => &
+            & l_continuity_mat(u_start(dim1):u_end(dim1),:)
+    end do
+
+    ! ( M    C  -L)(u)   (v)
+    ! ( -C^T N  0 )(h) = (j)
+    ! ( L^T  0  0 )(l)   (0)
+    ! 
+    ! (u)   (M    C)^{-1}(v)   (M    C)^{-1}(L)
+    ! (h) = (-C^T N)     (j) + (-C^T N)     (0)(l)
+    ! so
+    !        (M    C)^{-1}(L)         (M    C)^{-1}(v)
+    ! (L^T 0)(-C^T N)     (0)=-(L^T 0)(-C^T N)     (j)
+
+    !Get the local_solver matrix that obtains U and D from Lambda on the
+    !boundaries
+    call get_local_solver(local_solver_matrix,U,X,down,D,f,ele,&
+         & g,dt,theta,D0,pullback=.false.,weight=1.0,&
+         & have_constraint=have_constraint,&
+         & projection=.false.,poisson=.false.,&
+         & local_solver_rhs=local_solver_rhs)
+
+    !!!Construct the continuity matrix that multiplies lambda in 
+    !!! the U equation
+    !allocate l_continuity_mat
+    l_continuity_mat = 0.
+    !get list of neighbours
+    neigh => ele_neigh(D,ele)
+    !calculate l_continuity_mat
+    do ni = 1, size(neigh)
+       face=ele_face(U, ele, neigh(ni))
+       allocate(continuity_face_mat(mdim,face_loc(U,face)&
+            &,face_loc(lambda_rhs,face)))
+       continuity_face_mat = 0.
+       call get_continuity_face_mat(continuity_face_mat,face,&
+            U,lambda_rhs)
+       do dim1 = 1, mdim
+          continuity_mat_u_ptr(dim1)%ptr(face_local_nodes(U,face),&
+               face_local_nodes(lambda_rhs,face))=&
+          continuity_mat_u_ptr(dim1)%ptr(face_local_nodes(U,face),&
+               face_local_nodes(lambda_rhs,face))+&
+               continuity_face_mat(dim1,:,:)
+       end do
+       do  dim1 = 1, mdim
+          call addto(continuity_mat,dim1,1,face_global_nodes(U,face)&
+               &,face_global_nodes(lambda_rhs,face),&
+               &continuity_face_mat(dim1,:,:))
+       end do
+
+       deallocate(continuity_face_mat)
+    end do
+
+    !compute l_continuity_mat2 = inverse(local_solver)*l_continuity_mat
+    allocate(l_continuity_mat2(uloc*2+dloc+n_constraints,lloc))
+    l_continuity_mat2 = l_continuity_mat
+    call solve(local_solver_matrix,l_continuity_mat2)
+
+    !compute helmholtz_loc_mat
+    allocate(helmholtz_loc_mat(lloc,lloc))
+    helmholtz_loc_mat = matmul(transpose(l_continuity_mat),l_continuity_mat2)
+
+    !insert helmholtz_loc_mat into global lambda matrix
+    call addto(lambda_mat,ele_nodes(lambda_rhs,ele),&
+         ele_nodes(lambda_rhs,ele),helmholtz_loc_mat)
+  end subroutine assemble_newton_solver_ele
  
   subroutine assemble_hybridized_helmholtz_ele(D,f,U,X,down,ele, &
-       g,dt,theta,D0,lambda_mat,lambda_rhs,U_rhs,D_rhs,&
+       g,dt,theta,D0,pullback,weight,lambda_mat,lambda_rhs,U_rhs,D_rhs,&
        continuity_mat,projection,poisson,u_rhs_local)
     !subroutine to assemble hybridized helmholtz equation.
     !For assembly, must provide:
@@ -243,7 +705,8 @@
     type(vector_field), intent(in), optional :: U_rhs
     type(scalar_field), intent(in), optional :: D_rhs
     integer, intent(in) :: ele
-    real, intent(in) :: g,dt,theta,D0
+    logical, intent(in) :: pullback
+    real, intent(in) :: g,dt,theta,D0,weight
     type(csr_matrix), intent(inout) :: lambda_mat
     type(block_csr_matrix), intent(inout), optional :: continuity_mat
     logical, intent(in), optional :: projection, poisson,u_rhs_local
@@ -259,7 +722,7 @@
     real, dimension(:),allocatable,target :: Rhs_loc
     real, dimension(:,:), allocatable :: local_solver_matrix, local_solver_rhs
     type(element_type) :: U_shape
-    integer :: stat, d_start, d_end, dim1, mdim, uloc,dloc, lloc
+    integer :: stat, d_start, d_end, dim1, mdim, uloc,dloc, lloc, iloc
     integer, dimension(mesh_dim(U)) :: U_start, U_end
     type(real_vector), dimension(mesh_dim(U)) :: rhs_u_ptr
     real, dimension(:), pointer :: rhs_d_ptr
@@ -321,8 +784,15 @@
     !Get the local_solver matrix that obtains U and D from Lambda on the
     !boundaries
     call get_local_solver(local_solver_matrix,U,X,down,D,f,ele,&
-         & g,dt,theta,D0,have_constraint,local_solver_rhs,&
-         projection)
+         & g,dt,theta,D0,pullback,weight,have_constraint,&
+         & projection=projection,poisson=poisson,&
+         & local_solver_rhs=local_solver_rhs)
+    if(poisson.and.(ele==1)) then
+          !Fix the zero level of D
+       local_solver_matrix(d_start,:) = 0.
+       local_solver_matrix(:,d_start) = 0.
+       local_solver_matrix(d_start,d_start) = 1.
+    end if
 
     !!!Construct the continuity matrix that multiplies lambda in 
     !!! the U equation
@@ -343,13 +813,13 @@
                face_local_nodes(lambda_rhs,face))=&
           continuity_mat_u_ptr(dim1)%ptr(face_local_nodes(U,face),&
                face_local_nodes(lambda_rhs,face))+&
-               continuity_face_mat(dim1,:,:)
+               weight*continuity_face_mat(dim1,:,:)
        end do
        if(present(continuity_mat)) then
           do  dim1 = 1, mdim
              call addto(continuity_mat,dim1,1,face_global_nodes(U,face)&
                   &,face_global_nodes(lambda_rhs,face),&
-                  &continuity_face_mat(dim1,:,:))
+                  &weight*continuity_face_mat(dim1,:,:))
           end do
        end if
 
@@ -368,7 +838,12 @@
     !construct lambda_rhs
     rhs_loc=0.
     lambda_rhs_loc = 0.
-    call assemble_rhs_ele(Rhs_loc,D,U,X,ele,D_rhs,U_rhs,u_rhs_local)
+    call assemble_rhs_ele(Rhs_loc,D,U,X,ele,pullback,&
+         weight,D_rhs,U_rhs,u_rhs_local)
+    if(poisson.and.(ele==1)) then
+       !Fix the zero level of D
+       rhs_loc(d_start)=0.0
+    end if
     if(.not.(present(d_rhs).or.present(u_rhs)))then
        assert(.not.present_and_true(projection))
        assert(.not.present_and_true(poisson))
@@ -384,9 +859,9 @@
          ele_nodes(lambda_rhs,ele),helmholtz_loc_mat)
 
   end subroutine assemble_hybridized_helmholtz_ele
-
+  
   subroutine reconstruct_U_d_ele(D,f,U,X,down,ele, &
-       g,dt,theta,D0,U_rhs,D_rhs,lambda,&
+       g,dt,theta,D0,pullback,weight,U_rhs,D_rhs,lambda,&
        &D_out,U_out,projection,poisson,u_rhs_local)
     !subroutine to reconstruct U and D having solved for lambda
     implicit none
@@ -399,12 +874,13 @@
     type(scalar_field), intent(inout), optional :: D_out
     type(vector_field), intent(inout), optional :: U_out
     integer, intent(in) :: ele
-    real, intent(in) :: g,dt,theta,D0
+    logical, intent(in) :: pullback
+    real, intent(in) :: g,dt,theta,D0,weight
     logical, intent(in), optional :: projection, poisson,u_rhs_local
     !
     real, allocatable, dimension(:,:,:) :: continuity_face_mat
     integer :: ni, face
-    integer, dimension(:), pointer :: neigh
+    integer, dimension(:), pointer :: neigh, d_nodes
     type(element_type) :: U_shape
     integer :: d_start, d_end, dim1, mdim, uloc,dloc,lloc
     integer, dimension(mesh_dim(U)) :: U_start, U_end
@@ -420,7 +896,7 @@
     logical :: have_constraint
     integer :: n_constraints, i1
     type(constraints_type), pointer :: constraints
-    real :: constraint_check
+    real :: constraint_check, u_max
 
     !Get some sizes
     lloc = ele_loc(lambda,ele)
@@ -460,12 +936,23 @@
     !Get the local_solver matrix that obtains U and D from Lambda on the
     !boundaries
     call get_local_solver(local_solver_matrix,U,X,down,D,f,ele,&
-         & g,dt,theta,D0,have_constraint,&
-         & local_solver_rhs=local_solver_rhs,projection=projection,&
-         & poisson=poisson)
+         & g,dt,theta,D0,pullback,weight,have_constraint,&
+         & projection=projection,poisson=poisson,&
+         & local_solver_rhs=local_solver_rhs)
+    if(poisson.and.(ele==1)) then
+          !Fix the zero level of D
+       local_solver_matrix(d_start,:) = 0.
+       local_solver_matrix(:,d_start) = 0.
+       local_solver_matrix(d_start,d_start) = 1.
+    end if
 
     !Construct the rhs sources for U from lambda
-    call assemble_rhs_ele(Rhs_loc,D,U,X,ele,D_rhs,U_rhs,U_rhs_local)
+    call assemble_rhs_ele(Rhs_loc,D,U,X,ele,pullback,&
+         weight,D_rhs,U_rhs,U_rhs_local)
+    if(poisson.and.(ele==1)) then
+       !Fix the zero level of D
+       rhs_loc(d_start)=0.0
+    end if
     if(.not.(present(d_rhs).or.present(u_rhs)))then
        assert(.not.present_and_true(projection))
        assert(.not.present_and_true(poisson))
@@ -486,7 +973,7 @@
        do dim1 = 1, mdim
           rhs_u_ptr(dim1)%ptr(face_local_nodes(U,face)) = &
                & rhs_u_ptr(dim1)%ptr(face_local_nodes(U,face)) + &
-               & matmul(continuity_face_mat(dim1,:,:),&
+               & weight*matmul(continuity_face_mat(dim1,:,:),&
                &        face_val(lambda,face))
        end do
        deallocate(continuity_face_mat)
@@ -500,9 +987,6 @@
     ! (h) = (-C^T N)     (j) + (-C^T N)     (0)(l)
 
     call solve(local_solver_matrix,Rhs_loc)
-    rhs_loc = matmul(local_solver_matrix,rhs_loc)
-    call solve(local_solver_matrix,Rhs_loc)
-
     do dim1 = 1, mdim
        U_solved(dim1,:) = rhs_loc(u_start(dim1):u_end(dim1))
        if(.not.(present_and_true(poisson))) then
@@ -513,7 +997,7 @@
           end if
        end if
     end do
-    
+
     D_solved = rhs_loc(d_start:d_end)
     if(.not.(present_and_true(projection))) then
        if(present(D_out)) then
@@ -524,26 +1008,29 @@
     end if
 
     !check that the constraints are satisfied
+    !this is just a rescaling so don't need to use weights
     if(have_constraint) then
        constraints => U%mesh%shape%constraints
        do i1 = 1, constraints%n_constraints
           constraint_check = 0.
+          u_max = 0.
           do dim1 = 1, mdim
+             u_max = max(u_max,maxval(abs(U_solved(dim1,:))))
              constraint_check = constraint_check + &
                   & sum(U_solved(dim1,:)*constraints%orthogonal(i1,:,dim1))
           end do
-          if(abs(constraint_check)>1.0e-8) then
-             ewrite(1,*) 'Constraint check', constraint_check
+          if(abs(constraint_check)/(max(1.0,u_max))>1.0e-8) then
+             ewrite(2,*) 'Constraint check', constraint_check
              FLAbort('Constraint not enforced')
           end if
-       end do       
+       end do
     end if
 
   end subroutine reconstruct_U_d_ele
 
   subroutine get_local_solver(local_solver_matrix,U,X,down,D,f,ele,&
-       & g,dt,theta,D0,have_constraint, &
-       & local_solver_rhs,projection,poisson)
+       & g,dt,theta,D0,pullback,weight,have_constraint, &
+       & projection,poisson,local_solver_rhs)
     !Subroutine to get the matrix and rhs for obtaining U and D within
     !element ele from the lagrange multipliers on the boundaries.
     !This matrix-vector system is referred to as the "local solver" in the 
@@ -554,15 +1041,16 @@
     implicit none
     !If projection is present and true, set dt to zero and just project U 
     !into div-conforming space
-    real, intent(in) :: g,dt,theta,D0
+    real, intent(in) :: g,dt,theta,D0,weight
     type(vector_field), intent(in) :: U,X,down
     type(scalar_field), intent(in) :: D,f
     integer, intent(in) :: ele
+    logical, intent(in) :: pullback
     real, dimension(:,:)&
          &, intent(inout) :: local_solver_matrix
     real, dimension(:,:)&
          &, intent(inout), optional :: local_solver_rhs
-    logical, intent(in), optional :: projection, poisson
+    logical, intent(in) :: projection, poisson
     logical, intent(in) :: have_constraint
     !
     real, dimension(mesh_dim(U), X%dim, ele_ngi(U,ele)) :: J
@@ -583,12 +1071,20 @@
     integer, dimension(mesh_dim(U)) :: U_start, U_end
     type(constraints_type), pointer :: constraints
     integer :: i1, c_start, c_end
-    real :: l_dt
+    real :: l_dt,l_theta,l_d0,detJ_bar
 
-    if(present_and_true(projection)) then
+    if(projection) then
        l_dt = 0.
+       l_theta = 0.
+       l_d0 = 1.
+    else if(poisson) then
+       l_dt = 1.
+       l_theta = 1.
+       l_d0 = 1.
     else
        l_dt = dt
+       l_theta = theta
+       l_d0 = d0
     end if
 
     mdim = mesh_dim(U)
@@ -606,13 +1102,13 @@
     if(present(local_solver_rhs)) then
        local_solver_rhs = 0.
     end if
-    
+
     u_shape=ele_shape(u, ele)
     D_shape=ele_shape(d, ele)
     D_ele => ele_nodes(D, ele)
     U_ele => ele_nodes(U, ele)
-    
-    if(present_and_true(projection)) then
+
+    if(projection.or.poisson) then
        f_gi = 0.
     else
        f_gi = ele_val_at_quad(f,ele)
@@ -625,6 +1121,7 @@
     !detwei is needed for pressure mass matrix
     call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), J=J, &
          detwei=detwei, detJ=detJ)
+    detJ_bar = sum(detJ)/size(detJ)
 
     !----construct local solver
     !metrics for velocity mass and coriolis matrices
@@ -643,51 +1140,71 @@
     !dt*theta*<\phi,div u>         + D_0<\phi,h>         = 0 
 
     !pressure mass matrix (done in global coordinates)
-    local_solver_matrix(d_start:d_end,d_start:d_end)=&
-         &shape_shape(d_shape,d_shape,detwei)
+    !not included in pressure solver
+    if(.not.poisson) then
+       if(pullback) then
+          local_solver_matrix(d_start:d_end,d_start:d_end)=&
+               &shape_shape(d_shape,d_shape,detJ_bar**2/detJ*&
+               D_shape%quadrature%weight)
+       else
+          local_solver_matrix(d_start:d_end,d_start:d_end)=&
+               &shape_shape(d_shape,d_shape,detwei)
+       end if
        if(present(local_solver_rhs)) then
-          local_solver_rhs(d_start:d_end,d_start:d_end) = &
-               shape_shape(d_shape,d_shape,detwei)
+          if(pullback) then
+             local_solver_rhs(d_start:d_end,d_start:d_end) = &
+                  shape_shape(d_shape,d_shape,detJ_bar**2/detJ*&
+                  D_shape%quadrature%weight)
+          else
+             local_solver_rhs(d_start:d_end,d_start:d_end) = &
+                  shape_shape(d_shape,d_shape,detwei)
+          end if
        end if
+    end if
     !divergence matrix (done in local coordinates)
-    l_div_mat = dshape_shape(u_shape%dn,d_shape,&
-         &D_shape%quadrature%weight)
+    if(pullback) then
+       l_div_mat = dshape_shape(u_shape%dn,d_shape,&
+            &detJ_bar/detJ*D_shape%quadrature%weight)
+    else
+       l_div_mat = dshape_shape(u_shape%dn,d_shape,D_shape%quadrature%weight)
+    end if
     do dim1 = 1, mdim
        !pressure gradient term [integrated by parts so minus sign]
        local_solver_matrix(u_start(dim1):u_end(dim1),d_start:d_end)=&
-            & -g*l_dt*theta*l_div_mat(dim1,:,:)
+            & -g*l_dt*l_theta*weight*l_div_mat(dim1,:,:)
        if(present(local_solver_rhs)) then
           local_solver_rhs(u_start(dim1):u_end(dim1),d_start:d_end)=&
-               & -g*(theta-1.0)*l_dt*l_div_mat(dim1,:,:)
+               & -g*(l_theta-1.0)*l_dt*weight*l_div_mat(dim1,:,:)
        end if
        !divergence continuity term
        local_solver_matrix(d_start:d_end,u_start(dim1):u_end(dim1))=&
-            & d0*l_dt*theta*transpose(l_div_mat(dim1,:,:))
+            & l_d0*l_dt*l_theta*weight*transpose(l_div_mat(dim1,:,:))
        if(present(local_solver_rhs)) then
           local_solver_rhs(d_start:d_end,u_start(dim1):u_end(dim1))=&
-               & d0*(theta-1.0)*l_dt*transpose(l_div_mat(dim1,:,:))
+               & l_d0*(l_theta-1.0)*l_dt*weight*transpose(l_div_mat(dim1,:,:))
        end if
     end do
     !velocity mass matrix and Coriolis matrix (done in local coordinates)
     l_u_mat = shape_shape_tensor(u_shape, u_shape, &
-         u_shape%quadrature%weight, Metric+l_dt*theta*Metricf)
+         u_shape%quadrature%weight, Metric+l_dt*l_theta*Metricf)
+
     do dim1 = 1, mdim
        do dim2 = 1, mdim
           local_solver_matrix(u_start(dim1):u_end(dim1),&
                u_start(dim2):u_end(dim2))=&
-               & l_u_mat(dim1,dim2,:,:)
+               & weight**2*l_u_mat(dim1,dim2,:,:)
        end do
     end do
 
     if(present(local_solver_rhs)) then
        l_u_mat = shape_shape_tensor(u_shape, u_shape, &
             u_shape%quadrature%weight, &
-            Metric+(theta-1.0)*l_dt*Metricf)
+            Metric+(l_theta-1.0)*l_dt*Metricf)
        do dim1 = 1, mdim
           do dim2 = 1, mdim
              local_solver_rhs(u_start(dim1):u_end(dim1),&
                   u_start(dim2):u_end(dim2))=&
-                  & l_u_mat(dim1,dim2,:,:)
+                  & weight**2*l_u_mat(dim1,dim2,:,:)
           end do
        end do
     end if
@@ -708,52 +1225,6 @@
         
   end subroutine get_local_solver
 
-  subroutine get_up_gi(X,ele,up_gi,orientation)
-    !subroutine to replace up_gi with a normal to the surface
-    !with the same orientation
-    implicit none
-    type(vector_field), intent(in) :: X
-    integer, intent(in) :: ele
-    real, dimension(X%dim,ele_ngi(X,ele)), intent(inout) :: up_gi
-    integer, intent(out), optional :: orientation
-    !
-    real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J
-    integer :: gi
-    real, dimension(X%dim,ele_ngi(X,ele)) :: normal_gi
-    real, dimension(ele_ngi(X,ele)) :: orientation_gi
-    real :: norm
-
-    call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), J=J)
-
-    select case(mesh_dim(X)) 
-    case (2)
-       !Coriolis only makes sense for 2d surfaces embedded in 3d
-       do gi = 1, ele_ngi(X,ele)
-          normal_gi(:,gi) = cross_product(J(1,:,gi),J(2,:,gi))
-          norm = sqrt(sum(normal_gi(:,gi)**2))
-          normal_gi(:,gi) = normal_gi(:,gi)/norm
-       end do
-       do gi = 1, ele_ngi(X,ele)
-          orientation_gi(gi) = dot_product(normal_gi(:,gi),up_gi(:,gi))
-       end do
-       if(any(abs(orientation_gi-orientation_gi(1))>1.0e-8)) then
-          FLAbort('Nasty geometry problem')
-       end if
-       do gi = 1, ele_ngi(X,ele)
-          up_gi(:,gi) = normal_gi(:,gi)*orientation_gi(gi)
-       end do
-       if(present(orientation)) then
-          if(orientation_gi(1)>0.0) then
-             orientation = 1
-          else
-             orientation = -1
-          end if
-       end if
-    case default
-       FLAbort('not implemented')
-    end select
-  end subroutine get_up_gi
-
   function get_orientation(X_val, up) result (orientation)
     !function compares the orientation of the element with the 
     !up direction
@@ -853,73 +1324,6 @@
 
   end subroutine get_scalar_continuity_face_mat
 
-  subroutine get_local_normal(norm,weight,U,face)
-    !Function returns normal to face on local 2D element
-    implicit none
-    type(vector_field), intent(in) :: U
-    integer, intent(in) :: face
-    real, dimension(U%dim, face_ngi(U,face)), intent(out) :: norm
-    real, intent(out) :: weight
-
-    integer :: i
-
-    select case(U%mesh%shape%numbering%family)
-    case (FAMILY_SIMPLEX)
-       if(U%dim==1) then
-          if(face==1) then
-             forall(i=1:face_ngi(U,face)) norm(1,i)=1.
-          else if(face==2) then
-             forall(i=1:face_ngi(U,face)) norm(1,i)=-1.
-          else 
-             FLAbort('Funny face?')
-          end if
-          weight = 1.0
-
-       else if(U%dim==2) then
-          if(face==1) then
-             forall(i=1:face_ngi(U,face)) norm(1:2,i)=(/-1.,0./)
-          else if(face==2) then
-             forall(i=1:face_ngi(U,face)) norm(1:2,i)=(/0.,-1./)
-          else if(face==3) then
-             forall(i=1:face_ngi(U,face)) norm(1:2,i)=(/1/sqrt(2.),1&
-                  &/sqrt(2.)/)
-          else 
-             FLAbort('Funny face?')
-          end if
-
-          !Integral is taken on one of the edges of the local 2D element
-          !This edge must be transformed to the local 1D element
-          !to do numerical integration, with the following weight factors
-          if(face==3) then
-             weight = sqrt(2.)
-          else
-             weight = 1.0
-          end if
-
-       else
-          FLAbort('Dimension not supported.')
-       end if
-    case (FAMILY_CUBE)
-       if(U%dim==2) then
-          if(face==1) then
-             forall(i=1:face_ngi(U,face)) norm(1:2,i)=(/-1.,0./)
-          else if(face==2) then
-             forall(i=1:face_ngi(U,face)) norm(1:2,i)=(/ 1.,0./)
-          else if(face==3) then
-             forall(i=1:face_ngi(U,face)) norm(1:2,i)=(/0.,-1./)
-          else if(face==4) then
-             forall(i=1:face_ngi(U,face)) norm(1:2,i)=(/0.,1./)
-          else
-             FLAbort('Funny face?')
-          end if
-          weight = 1.0
-       else
-          FLAbort('Dimension not supported.')
-       end if
-    end select
-
-  end subroutine get_local_normal
-
   subroutine compute_cartesian_ele(U_cart,U,X,ele)
     implicit none
     type(vector_field), intent(inout) :: U_cart
@@ -935,20 +1339,20 @@
     integer, dimension(mesh_dim(U)) :: U_start, U_end
     integer, dimension(:), pointer :: U_ele
     integer :: mdim, uloc, gi
-    real, dimension(ele_ngi(U,ele)) :: detwei, detJ
+    real, dimension(ele_ngi(U,ele)) :: detwei, detJ 
     real, dimension(mesh_dim(U), X%dim, ele_ngi(U,ele)) :: J
 
-    mdim = mesh_dim(U)
-    uloc = ele_loc(U,ele)
-    do dim1 = 1, mdim
-       u_start(dim1) = uloc*(dim1-1)+1
-       u_end(dim1) = uloc+dim1
+    mdim = mesh_dim(U) 
+    uloc = ele_loc(U,ele) 
+    do dim1 = 1, mdim 
+       u_start(dim1) = uloc*(dim1-1)+1 
+       u_end(dim1) = uloc+dim1 
     end do
 
     U_ele => ele_nodes(U, ele)
 
-    u_shape=ele_shape(u, ele)
-    call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), J=J, &
+    u_shape=ele_shape(u, ele) 
+    call compute_jacobian(ele_val(X,ele),ele_shape(X,ele), J=J, & 
          detwei=detwei,detJ=detJ)
 
     local_u_gi = ele_val_at_quad(U,ele)
@@ -1002,15 +1406,16 @@
     u1 = face_val_at_quad(U,face)
     u2 = face_val_at_quad(U,face2)
     jump = maxval(abs(sum(u1*n1+u2*n2,1)))
-    ewrite(1,*) jump
+    ewrite(2,*) jump
     assert(jump<1.0e-8)
 
   end subroutine check_continuity_local_face
 
-  subroutine check_continuity_ele(U_cart,X,ele)
+  subroutine check_continuity_ele(U_cart,X,ele,tolerance)
     implicit none
     type(vector_field), intent(in) :: U_cart,X
     integer, intent(in) :: ele
+    real, intent(in) :: tolerance
     !
     integer, dimension(:), pointer :: neigh
     integer :: ni,face,ele2,face2
@@ -1024,18 +1429,20 @@
        else
           face2 = -1
        end if
-       call check_continuity_face(U_cart,X,ele,ele2,face,face2)
+       call check_continuity_face(U_cart,X,ele,ele2,face,face2,tolerance)
     end do
   end subroutine check_continuity_ele
 
-  subroutine check_continuity_face(U_cart,X,ele,ele2,face,face2)
+  subroutine check_continuity_face(U_cart,X,ele,ele2,face,face2,tolerance)
     !subroutine to check the continuity of normal component
     !of velocity at quadrature points
     implicit none
     type(vector_field), intent(in) :: U_cart,X
     integer, intent(in) :: face,face2,ele,ele2
+    real, intent(in) :: tolerance
+    !
     real, dimension(X%dim, face_ngi(U_cart, face)) :: n1,n2
-    real, dimension(X%dim, face_ngi(U_cart, face)) :: u1,u2,x1
+    real, dimension(X%dim, face_ngi(U_cart, face)) :: u1,u2,x1,x2
     real, dimension(face_ngi(U_cart, face)) :: jump_at_quad
     integer :: dim1
     !
@@ -1043,6 +1450,13 @@
     x1 = face_val_at_quad(X,face)
     if(ele2>0) then
        u2 = face_val_at_quad(U_cart,face2)
+       x2 = face_val_at_quad(X,face)
+       if(any(x1.ne.x2)) then
+          ewrite(0,*) 'Face 1 X', x1
+          ewrite(0,*) 'Face 2 X', x2
+          FLExit('Something wrong with mesh?')
+       end if
+
     else
        u2 = 0.
     end if
@@ -1054,98 +1468,22 @@
        n2 = -n1
     end if
     jump_at_quad = sum(n1*u1+n2*u2,1)
-    if(maxval(abs(jump_at_quad))>1.0e-8) then
-       ewrite(1,*) 'Jump at quadrature face, face2 =', jump_at_quad
-       ewrite(1,*) 'ELE = ',ele,ele2
+    if(maxval(abs(jump_at_quad))>tolerance) then
+       ewrite(2,*) 'Jump at quadrature face, face2 =', jump_at_quad
+       ewrite(2,*) 'ELE = ',ele,ele2
        do dim1 = 1, X%dim
-          ewrite(1,*) 'normal',dim1,n1(dim1,:)
-          ewrite(1,*) 'normal',dim1,n2(dim1,:)
-          ewrite(1,*) 'X',dim1,x1(dim1,:)
+          ewrite(2,*) 'normal',dim1,n1(dim1,:)
+          ewrite(2,*) 'normal',dim1,n2(dim1,:)
+          ewrite(2,*) 'X',dim1,x1(dim1,:)
        end do
-       ewrite(1,*) 'n cpt1',sum(n1*u1,1)
-       ewrite(1,*) 'n cpt2',sum(n2*u2,1)
-       ewrite(1,*) jump_at_quad/max(maxval(abs(u1)),maxval(abs(u2)))
+       ewrite(2,*) 'n cpt1',sum(n1*u1,1)
+       ewrite(2,*) 'n cpt2',sum(n2*u2,1)
+       ewrite(2,*) jump_at_quad/max(maxval(abs(u1)),maxval(abs(u2)))
        FLAbort('stopping because of jumps')
     end if
     
   end subroutine check_continuity_face
 
-  function get_face_normal_manifold(X,ele,face) result (normal)
-    implicit none
-    type(vector_field), intent(in) :: X
-    integer, intent(in) :: ele, face
-    real, dimension(X%dim,face_ngi(X,face)) :: normal
-    !
-    real, dimension(ele_ngi(X,ele)) :: detwei
-    real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J
-    real, dimension(face_ngi(X,face)) :: detwei_f
-    real, dimension(mesh_dim(X)-1, X%dim, face_ngi(X,face)) :: J_f
-    real, dimension(ele_loc(X,ele),ele_loc(X,ele)) :: X_mass_mat
-    real, dimension(mesh_dim(X), X%dim, ele_loc(X,ele)) :: J_loc
-    real, dimension(ele_loc(X,ele)) :: J_loc_rhs
-    real, dimension(mesh_dim(X), X%dim, face_ngi(X,face)) :: J_face_gi
-    integer :: dim1, dim2, gi
-    real, dimension(X%dim) :: ele_normal_gi, edge_tangent_gi, X_mid_ele,&
-         &X_mid_face
-    type(element_type) :: X_shape, X_face_shape
-    real, dimension(X%dim,ele_loc(X,ele)) :: X_ele
-    real, dimension(X%dim,face_loc(X,face)) :: X_face
-
-    call compute_jacobian(ele_val(X,ele),ele_shape(X,ele), J=J, &
-         detwei=detwei)
-    call compute_jacobian(face_val(X,face),face_Shape(X,face), J=J_f, &
-         detwei=detwei_f)
-
-    !Jacobian can be expanded without error in X function space
-    !so we map it to the basis function DOFs by projection
-    X_shape = ele_shape(X,ele)
-    X_face_shape = face_shape(X,face)
-    X_mass_mat = shape_shape(X_shape,X_shape,detwei)
-    do dim1 = 1, mesh_dim(X)
-       do dim2 = 1, X%dim
-          J_loc_rhs = shape_rhs(X_shape,J(dim1,dim2,:)*detwei)
-          call solve(X_mass_mat,J_loc_rhs)
-          J_loc(dim1,dim2,:) = J_loc_rhs
-       end do
-    end do
-
-    do dim1 = 1, mesh_dim(X)
-       do dim2 = 1, X%dim
-          J_face_gi(dim1,dim2,:) = &
-               & matmul(transpose(X_face_shape%n),&
-               J_loc(dim1,dim2,face_local_nodes(X,face)))
-       end do
-    end do
-
-    X_mid_ele = sum(ele_val(X,ele),2)/size(ele_val(X,ele),2)
-    X_mid_face = sum(face_val(X,face),2)/size(face_val(X,face),2)
-    select case(X%dim)
-    case (3)
-       select case (mesh_dim(X))
-       case (2)
-          do gi = 1, face_ngi(X,face)
-             !Get normal to element e on face quad points
-             ele_normal_gi = cross_product(J_face_gi(1,:,gi),&
-                  &J_face_gi(2,:,gi))
-             ele_normal_gi = ele_normal_gi/(norm2(ele_normal_gi))
-             !Get tangent to face f
-             edge_tangent_gi = J_f(1,:,gi)
-             edge_tangent_gi = edge_tangent_gi/norm2(edge_tangent_gi)
-             !Compute normal to face f in manifold
-             normal(:,gi) = cross_product(ele_normal_gi,edge_tangent_gi)
-             if(dot_product(normal(:,gi),X_mid_face-X_mid_ele)<0) then
-                normal(:,gi)=-normal(:,gi)
-             end if
-          end do
-       case default
-          FLAbort('dimension combination not implemented')
-       end select
-    case default
-       FLAbort('dimension combination not implemented')
-    end select
-
-  end function get_face_normal_manifold
-
   subroutine reconstruct_lambda_nc(lambda,lambda_nc,X,ele)
     type(scalar_field), intent(in) :: lambda
     type(scalar_field), intent(inout) :: lambda_nc
@@ -1203,9 +1541,12 @@
          &shape_shape(lambda_nc_face_shape,lambda_nc_face_shape,detwei)
   end subroutine get_nc_rhs_face
 
-  subroutine assemble_rhs_ele(Rhs_loc,D,U,X,ele,D_rhs,U_rhs,u_rhs_local)
+  subroutine assemble_rhs_ele(Rhs_loc,D,U,X,ele,pullback,&
+       weight,D_rhs,U_rhs,u_rhs_local)
     implicit none
     integer, intent(in) :: ele
+    real, intent(in) :: weight
+    logical, intent(in) :: pullback
     type(scalar_field), intent(in), optional, target :: D_rhs
     type(vector_field), intent(in), optional, target :: U_rhs
     type(vector_field), intent(in) :: X,U
@@ -1226,6 +1567,7 @@
     type(scalar_field), pointer :: l_d_rhs
     type(vector_field), pointer :: l_u_rhs
     real, dimension(mesh_dim(U), mesh_dim(U)) :: Metric
+    real :: detJ_bar
 
     !Get some sizes
     mdim = mesh_dim(U)
@@ -1245,6 +1587,7 @@
     
     if(.not.(present(D_rhs).or.present(u_rhs))) then
        !We are in timestepping mode.
+       !This will be multiplied by the local_solver_rhs matrix later.
        rhs_loc(d_start:d_end) = ele_val(D,ele)
        u_rhs_loc = ele_val(U,ele)
        do dim1 = 1, mdim
@@ -1259,27 +1602,38 @@
        
        call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), J=J, &
             detJ=detJ,detwei=detwei)
+       detJ_bar = sum(detJ)/size(detJ)
        
        Rhs_loc = 0.
        if(have_d_rhs) then
-          Rhs_loc(d_start:d_end) = shape_rhs(ele_shape(D,ele),&
-               &ele_val_at_quad(l_D_rhs,ele)*detwei)
+          if(pullback) then
+             Rhs_loc(d_start:d_end) = shape_rhs(ele_shape(D,ele),&
+                  &ele_val_at_quad(l_D_rhs,ele)*&
+                  &detJ_bar**2/detJ*u_shape%quadrature%weight)
+          else
+             Rhs_loc(d_start:d_end) = shape_rhs(ele_shape(D,ele),&
+                  &ele_val_at_quad(l_D_rhs,ele)*detwei)
+          end if
        end if
        if(have_u_rhs) then
           if(present_and_true(u_rhs_local)) then
+             !Using local U, multiply by weight**2
              u_local_quad = ele_val_at_quad(l_u_rhs,ele)
              do gi=1,ele_ngi(U,ele)
                 Metric=matmul(J(:,:,gi), transpose(J(:,:,gi)))&
                      &/detJ(gi)
-                u_local_quad(:,gi) = matmul(Metric,u_local_quad(:,gi))
+                u_local_quad(:,gi) = matmul(Metric,u_local_quad(:,gi))&
+                     &*weight**2
              end do
           else
+             !Using cartesian U, multiply by weight
              allocate(u_cart_quad(l_U_rhs%dim,ele_ngi(X,ele)))
              u_cart_quad = ele_val_at_quad(l_u_rhs,ele)
              do gi = 1, ele_ngi(D,ele)
-                !Don't divide by detJ as we can use weight instead of detwei
+                !Don't divide by detJ as we can use quadrature weight
+                !instead of detwei
                 u_local_quad(:,gi) = matmul(J(:,:,gi)&
-                     &,u_cart_quad(:,gi))
+                     &,u_cart_quad(:,gi))*weight
              end do
           end if
           U_rhs_loc = shape_vector_rhs(u_shape,&
@@ -1361,8 +1715,8 @@
        call compute_energy_ele(energy,U,D,X,D0,g,ele)
     end do
 
-    ewrite(1,*) 'Energy:= ', energy
-    ewrite(1,*) 'Change in energy:= ', energy-old_energy
+    ewrite(2,*) 'Energy:= ', energy
+    ewrite(2,*) 'Percentage Change in energy:= ', (energy-old_energy)/energy
 
   end subroutine compute_energy_hybridized
 
@@ -1414,12 +1768,14 @@
     type(state_type), intent(inout) :: state
     !
     type(scalar_field), pointer :: D,psi,f
-    type(scalar_field) :: D_rhs
+    type(scalar_field) :: D_rhs,tmp_field
     type(vector_field), pointer :: U_local,down,X, U_cart
-    type(vector_field) :: Coriolis_term, Balance_eqn, tmp_field
-    integer :: ele,dim1
-    real :: g
-    logical :: elliptic_method
+    type(vector_field) :: Coriolis_term, Balance_eqn, tmpV_field
+    integer :: ele,dim1,i1, stat
+    real :: g, D0
+    logical :: elliptic_method, pullback
+    real :: u_max, b_val, h_mean, area
+    real, dimension(:), allocatable :: weights
 
     D=>extract_scalar_field(state, "LayerThickness")
     psi=>extract_scalar_field(state, "Streamfunction")
@@ -1429,133 +1785,135 @@
     X=>extract_vector_field(state, "Coordinate")
     down=>extract_vector_field(state, "GravityDirection")
     call get_option("/physical_parameters/gravity/magnitude", g)
-    call allocate(tmp_field,mesh_dim(U_local), U_local%mesh, "tmp_field")
-       call allocate(D_rhs,D%mesh,'BalancedSolverRHS')
+    call allocate(tmpV_field,mesh_dim(U_local), U_local%mesh, "tmpV_field")
+    call allocate(D_rhs,D%mesh,'BalancedSolverRHS')
+    call allocate(tmp_field,D%mesh,'TmpField')
     call allocate(Coriolis_term,mesh_dim(U_local),&
          U_local%mesh,"CoriolisTerm")
     call allocate(balance_eqn,mesh_dim(D),u_local%mesh,'BalancedEquation')
 
+    !compute rescaling weights for local coordinates
+    allocate(weights(element_count(X)))
+    weights = 1.0
+    !call get_weights(X,weights)
+
     !STAGE 1: Set velocity from streamfunction
     do ele = 1, element_count(D)
        call set_local_velocity_from_streamfunction_ele(&
-            &U_local,psi,down,X,ele)
+            &U_local,psi,down,X,ele,weights(ele))
     end do
 
     !STAGE 1a: verify that velocity projects is div-conforming
-    call project_local_to_cartesian(X,U_local,U_cart)
-    do ele = 1, ele_count(U_local)
-       call check_continuity_ele(U_cart,X,ele)
-    end do
-
+    !call project_local_to_cartesian(X,U_local,U_cart,weights=weights)
+    !do ele = 1, ele_count(U_local)
+    !   call check_continuity_ele(U_cart,X,ele,tolerance=1.0e-8)
+    !end do
     !Stage 1b: verify that projection is idempotent
-    ewrite(1,*) 'CHECKING CONTINUOUS', maxval(abs(u_local%val))
-    call solve_hybridized_helmholtz(state,U_Rhs=U_local,&
-         &U_out=tmp_field,&
-         &compute_cartesian=.true.,&
-         &check_continuity=.true.,projection=.true.,&
-         &u_rhs_local=.true.)!verified that projection is idempotent
-    assert(maxval(abs(U_local%val-tmp_field%val))<1.0e-8)
-
-    elliptic_method = .false.
-
-    if(elliptic_method) then
-
-       !STAGE 2: Construct Coriolis term
+    ewrite(2,*) 'CHECKING CONTINUOUS', maxval(abs(u_local%val))
+    tmpV_field%val = U_local%val
+    call project_to_constrained_space(state,tmpV_field)
+    ewrite(2,*) maxval(abs(U_local%val-tmpV_field%val)), 'continuity'
+    u_max = 0.
+    do dim1 = 1, mesh_dim(D)
+       u_max = max(u_max,maxval(abs(U_local%val)))
+    end do
+    assert(maxval(abs(U_local%val-tmpV_field%val)/max(1.0,u_max))<1.0e-8)
+
+    if(have_option("/material_phase::Fluid/vector_field::Velocity/prognostic&
+         &/initial_condition::WholeMesh/balanced/elliptic_solver")) then
+       ewrite(2,*) 'Using elliptic solver'
+       !Construct Coriolis term
        call zero(Coriolis_term)
        do ele = 1, element_count(D)
+          !weights not needed since both sides of equation are multiplied
+          !by weight(ele)^2 in each element
           call set_coriolis_term_ele(Coriolis_term,f,down,U_local,X,ele)
-       end do!checked!signs checked
-
-       !STAGE 3: Project Coriolis term into div-conforming space
-
-       !debugging bits - checking if it works with cartesian instead
-       call project_local_to_cartesian(X,Coriolis_Term,U_cart)
-       call solve_hybridized_helmholtz(state,U_Rhs=U_cart,&
-            &U_out=tmp_field,&
-            &compute_cartesian=.true.,output_dense=.false.,&
-            &check_continuity=.true.,projection=.true.,&
-            &u_rhs_local=.false.)
-
-       ewrite(0,*) 'REMEMBER TO REMOVE DEBUGGING TESTS'
-       call solve_hybridized_helmholtz(state,U_Rhs=Coriolis_term,&
-            &U_out=tmp_field,&
-            &compute_cartesian=.true.,&
-            &check_continuity=.true.,projection=.true.,&
-            &u_rhs_local=.true.)!verified that projection is idempotent
-
-       !STAGE 4: Construct the RHS for the balanced layer depth equation
-       call zero(D_rhs)
-       do ele = 1, element_count(D)
-          call set_geostrophic_balance_rhs_ele(D_rhs,Coriolis_term,ele)
-       end do
-
-       !STAGE 5: Solve Poisson equation for the balanced layer depth
-       ewrite(0,*) 'REMEMBER ABOUT SETTING MEAN VALUE'
-       ewrite(0,*) trim(u_cart%option_path)
-
-       call solve_hybridized_helmholtz(state,D_rhs=D_rhs,&
-            &compute_cartesian=.false.,&
-            &check_continuity=.false.,Poisson=.true.,&
-            &solver_option_path=trim(u_cart%option_path)//'/prognostic/initial_condition::WholeMesh/balanced')
-
-       !STAGE 6: Check if we have a balanced solution
-       !Can be done by projecting balance equation into div-conforming space
-       !and checking that it is equal to zero
-       !STAGE 6a: Project balance equation into DG space
-       do ele = 1, element_count(D)
-          call set_pressure_force_ele(balance_eqn,D,X,g,ele)
-       end do
-       call addto(balance_eqn,coriolis_term)
-
-       !STAGE 6b: Project balance equation into div-conforming space
-       call solve_hybridized_helmholtz(state,U_Rhs=balance_eqn,&
-            &U_out=balance_eqn,&
-            &compute_cartesian=.true.,&
-            &check_continuity=.true.,projection=.true.,&
-            &u_rhs_local=.true.)
-
-       do dim1 = 1, mesh_dim(D)
-          ewrite(1,*) 'Balance equation', maxval(abs(balance_eqn%val(dim1,:)))
-          assert(maxval(abs(balance_eqn%val(dim1,:)))<1.0e-8)
-       end do
-
+       end do
+
+       !Project Coriolis term into div-conforming space
+       call project_to_constrained_space(state,Coriolis_term)
+
+       !Calculate divergence of Coriolis term 
+       call zero(d_rhs)
+
+       pullback = have_option('/material_phase::Fluid/scalar_field::LayerThickn&
+            &ess/prognostic/spatial_discretisation/discontinuous_galerkin/wave_&
+            &equation/pullback')
+       
+       do ele = 1, element_count(D)
+          call compute_divergence_ele(Coriolis_term,d_rhs,X,ele,pullback)
+       end do
+
+       !Solve for balanced layer depth
+
+       ewrite(2,*) 'Solving elliptic problem for balanced layer depth.'
+       call solve_hybridized_helmholtz(state,d_Rhs=d_rhs,&
+            &U_out=tmpV_field,d_out=tmp_field,&
+            &compute_cartesian=.true.,&
+            &poisson=.true.,u_rhs_local=.true.)
+       D%val = tmp_field%val
     else
+       ewrite(2,*) 'Using streamfunction projection'
        !Project the streamfunction into pressure space
        do ele = 1, element_count(D)
           call project_streamfunction_for_balance_ele(D,psi,X,f,g,ele)
        end do
-       ewrite(1,*) maxval(abs(D%val))
-
-       !debugging tests
-       call zero(Coriolis_term)
-       do ele = 1, element_count(D)
-          call set_coriolis_term_ele(Coriolis_term,f,down,U_local,X,ele)
-       end do
-       call zero(balance_eqn)
-       do ele = 1, element_count(D)
-          call set_pressure_force_ele(balance_eqn,D,X,g,ele)
-       end do
-       call addto(balance_eqn,coriolis_term,scale=1.0)
-       ewrite(1,*) 'CJC b4',maxval(abs(balance_eqn%val)),&
-            & maxval(abs(coriolis_term%val))
-       !Project balance equation into div-conforming space
-       call solve_hybridized_helmholtz(state,U_Rhs=balance_eqn,&
-            &U_out=balance_eqn,&
-            &compute_cartesian=.true.,&
-            &check_continuity=.true.,projection=.true.,&
-            &u_rhs_local=.true.)
-
-        do dim1 = 1, mesh_dim(D)
-           ewrite(1,*) 'Balance equation', maxval(abs(balance_eqn%val(dim1,:)))
-           assert(maxval(abs(balance_eqn%val(dim1,:)))<1.0e-8)
-        end do
     end if
-    !Clean up after yourself
+
+    !Subtract off the mean part
+    !Didn't bother with pullback
+    h_mean = 0.
+    area = 0.
+    do ele = 1, element_count(D)
+       call assemble_mean_ele(D,X,h_mean,area,ele)
+    end do
+    h_mean = h_mean/area
+    D%val = D%val - h_mean
+    !Add back on the correct mean depth
+    call get_option("/material_phase::Fluid/scalar_field::LayerThickness/&
+         &prognostic/mean_layer_thickness",D0)
+    D%val = D%val + D0
+
+    !debugging tests
+    call zero(Coriolis_term)
+    do ele = 1, element_count(D)
+       !weights not needed since both sides of equation are multiplied
+       !by weight(ele)^2 in each element
+       call set_coriolis_term_ele(Coriolis_term,f,down,U_local,X,ele)
+    end do
+ 
+    call zero(balance_eqn)
+    do ele = 1, element_count(D)
+       call set_pressure_force_ele(balance_eqn,D,X,g,ele,weights(ele),&
+            pullback)
+    end do
+    call addto(balance_eqn,coriolis_term,scale=1.0)
+    ewrite(2,*) 'Project balance equation into div-conforming space'
+    ewrite(2,*) 'CJC b4',maxval(abs(balance_eqn%val)),&
+         & maxval(abs(coriolis_term%val))
+    b_val = maxval(abs(balance_eqn%val))
+    !Project balance equation into div-conforming space
+    call solve_hybridized_helmholtz(state,U_Rhs=balance_eqn,&
+         &U_out=balance_eqn,&
+         &compute_cartesian=.true.,&
+         &projection=.true.,&
+         &poisson=.false.,&
+         &u_rhs_local=.true.)
+    
+    do dim1 = 1, mesh_dim(D)
+       ewrite(2,*) 'Balance equation', maxval(abs(balance_eqn%val(dim1,:)))/b_val
+       !assert(maxval(abs(balance_eqn%val(dim1,:)))/b_val<1.0e-8)
+    end do
+    
+ !Clean up after yourself
     call deallocate(Coriolis_term)
     call deallocate(D_rhs)
     call deallocate(balance_eqn)
+    call deallocate(tmpV_field)
     call deallocate(tmp_field)
 
+    deallocate(weights)
+
   end subroutine set_velocity_from_geostrophic_balance_hybridized
 
   subroutine project_streamfunction_for_balance_ele(D,psi,X,f,g,ele)
@@ -1587,19 +1945,18 @@
 
   end subroutine project_streamfunction_for_balance_ele
   
-  subroutine set_pressure_force_ele(force,D,X,g,ele)
+  subroutine set_pressure_force_ele(force,D,X,g,ele,weight,pullback)
     implicit none
     type(vector_field), intent(inout) :: force
     type(scalar_field), intent(in) :: D
     type(vector_field), intent(in) :: X
-    real, intent(in) :: g
+    real, intent(in) :: g, weight
     integer, intent(in) :: ele
+    logical, intent(in) :: pullback
     !
     real, dimension(ele_ngi(D,ele)) :: D_gi
     real, dimension(mesh_dim(D),ele_loc(force,ele)) :: &
          & rhs_loc
-    real, dimension(ele_loc(force,ele),ele_loc(force,ele)) :: &
-         & l_mass_mat    
     integer :: dim1,dim2,gi,uloc
     real, dimension(mesh_dim(force), X%dim, ele_ngi(force,ele)) :: J
     real, dimension(ele_ngi(force,ele)) :: detJ
@@ -1609,6 +1966,7 @@
          mesh_dim(force)*ele_loc(force,ele)) :: l_u_mat
     real, dimension(mesh_dim(force)*ele_loc(force,ele)) :: force_rhs
     type(element_type) :: force_shape
+    real :: detJ_bar
 
     uloc = ele_loc(force,ele)
     force_shape = ele_shape(force,ele)
@@ -1618,11 +1976,14 @@
        Metric(:,:,gi)=matmul(J(:,:,gi), transpose(J(:,:,gi)))/detJ(gi)
     end do
 
+    detJ_bar = sum(detJ)/size(detJ)
+
     D_gi = ele_val_at_quad(D,ele)
+    if(pullback) then
+       D_gi = D_gi*detJ_bar/detJ
+    end if
     rhs_loc = -g*dshape_rhs(force%mesh%shape%dn,&
          D_gi*D%mesh%shape%quadrature%weight)
-    l_mass_mat = shape_shape(ele_shape(force,ele),ele_shape(force,ele),&
-         &force%mesh%shape%quadrature%weight)
     do dim1 = 1, mesh_dim(force)
        force_rhs((dim1-1)*uloc+1:dim1*uloc) = rhs_loc(dim1,:)
     end do
@@ -1635,10 +1996,12 @@
        end do
     end do
 
+    !LHS is multiplied by weights(ele)^2, RHS is multiplied by weights(ele)
+    
     call solve(l_u_mat,force_rhs)
     do dim1= 1, mesh_dim(force)
        call set(force,dim1,ele_nodes(force,ele),&
-            &force_rhs((dim1-1)*uloc+1:dim1*uloc))
+            &force_rhs((dim1-1)*uloc+1:dim1*uloc)/weight)
     end do
 
   end subroutine set_pressure_force_ele
@@ -1709,7 +2072,6 @@
     up_gi = -ele_val_at_quad(down,ele)
 
     call get_up_gi(X,ele,up_gi)
-
     coriolis_rhs = 0.
     l_u_mat = 0.
     !metrics for velocity mass and coriolis matrices
@@ -1747,15 +2109,16 @@
   end subroutine set_coriolis_term_ele
   
   subroutine set_local_velocity_from_streamfunction_ele(&
-       &U_local,psi,down,X,ele)
+       &U_local,psi,down,X,ele,weight)
     implicit none
     type(vector_field), intent(inout) :: U_local
     type(vector_field), intent(in) :: down,X
     type(scalar_field), intent(in) :: psi
     integer, intent(in) :: ele
+    real, intent(in) :: weight
     !
     real, dimension(ele_loc(psi,ele)) :: psi_loc
-    real, dimension(mesh_dim(psi),ele_ngi(psi,ele)) :: dpsi_gi
+    real, dimension(mesh_dim(psi),ele_ngi(psi,ele)) :: dpsi_gi, Uloc_gi
     real, dimension(ele_ngi(psi,ele)) :: div_gi
     real, dimension(mesh_dim(U_local),ele_loc(U_local,ele)) :: U_loc
     real, dimension(ele_loc(U_local,ele),ele_loc(U_local,ele)) :: &
@@ -1764,6 +2127,8 @@
     integer :: dim1,gi,uloc
     real, dimension(X%dim, ele_ngi(X,ele)) :: up_gi
     integer :: orientation
+    real, dimension(mesh_dim(U_local), X%dim, ele_ngi(U_local,ele)) :: J
+    real, dimension(ele_ngi(X,ele)) :: detJ
 
     uloc = ele_loc(U_local,ele)