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

[Richard Ferrier]

The darcy_impes program solves multiphase flow through porous media using an implicit pressure, explicit saturation (IMPES) method. See, for example, Mostaghini P., Tollit B.S., Neethling S.J., Gorman G.J. and Pain C.C., A control volume finite element method for adaptive mesh simulation of flow in heap leaching, J Eng Math, 1--11, 2013.

The test/darcy_impes_p1_* folders making up the test suite contain a number of different types of test (contact me for the details - I tried to insert a table here but the formatting got messed up). Basic tests include checks that the solver has converged, field values do not exceed certain boundary conditions, and so on. Regression tests are comparisons with previously computed solutions. Analytical comparisons are the most reliable checks that the program is behaving as expected, with comparisons made against a 1D quasi-analytical profile.

Brendan Tollit (btollit) implemented an extensive coverage of the three types of test, to which I have made a few changes. The general idea is to promote analytical comparisons over regression tests (or implement the latter only where analytical solutions are not available), since the regression solutions are sensitive to e.g. mesh design and are not guaranteed to be correct. It is intended that analytical tests will eventually come to replace regression tests, although I have erred on the side of caution and not deleted regression tests at this stage.

[Stephan Kramer (s-kramer)]

Hello, only had a quick look at this merge. Most code changes to code shared with fluidity seem reasonable at first glance.

I have some worries about the need to add this to trunk Fluidity, in terms of maintainability, etc. but these are probably better discussed in another forum.

The main problem at the moment is that this merge seems to add about 150M to a clean checkout (current checkout size is 300M). The main cause of this seem to be the tests/Tarantula_mesh/ and tests/Tarantula_mesh_v2 directories that were added. In general there seems to be a whole lot of things in there (reports and stuff) that are not appropriate in the tests directory. If such a large data set is needed for testing, I suggest changing it to a long test. Also, are we sure that there are no license issues associated with any of it? If we do remove these directories, we have to rebuild the history after their commits to avoid polluting the history with these.

[Richard Ferrier (rjferrier)]

Good point Stephen. I had completely missed the Tarantula folders. I will have a look and see what can be done with them. It will probably need discussion with Peyman who is on leave at the moment.

[James Robert Percival (j-percival)]

It appears that neither directory actually contains a test (in the test harness sense) at the moment anyway. As a secondary part of the question of the trunk merge, I note that the darcy solver gets built as default by the root makefile. Other similar codes (e.g. Colin Cotter's shallow water code) are only built if requested, or if the code is being tested.

[Gerard (g-gorman)]

Hi Richards

Can I suggest that you create a new branch to cache this tarantula work (call it tarantula), and then delete it entirely from the Darcy branch.

Cheers
Gerard

On 24 Feb 2014, at 15:06, Stephan Kramer <email address hidden> wrote:

> Hello, only had a quick look at this merge. Most code changes to code shared with fluidity seem reasonable at first glance.
>
> I have some worries about the need to add this to trunk Fluidity, in terms of maintainability, etc. but these are probably better discussed in another forum.
>
> The main problem at the moment is that this merge seems to add about 150M to a clean checkout (current checkout size is 300M). The main cause of this seem to be the tests/Tarantula_mesh/ and tests/Tarantula_mesh_v2 directories that were added. In general there seems to be a whole lot of things in there (reports and stuff) that are not appropriate in the tests directory. If such a large data set is needed for testing, I suggest changing it to a long test. Also, are we sure that there are no license issues associated with any of it? If we do remove these directories, we have to rebuild the history after their commits to avoid polluting the history with these.
> --
> https://code.launchpad.net/~fluidity-core/fluidity/darcy_impes_proper/+merge/207924
> Your team Fluidity Core Team is subscribed to branch lp:fluidity.

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

4136. By Richard Ferrier on 2014-02-24

Default building of darcy_impes binary removed from the root Makefile. The 3phase test folder's meshes had not been checked into R4133 (oops); they are herewith restored. This will be the last revision to contain the Tarantula test folders, deemed to bulky/experimental to include in the proposed merge. A new development branch will be created to keep those folders.

4137. By Richard Ferrier on 2014-02-24

Tarantula folders deleted.

4138. By Richard Ferrier on 2014-02-25

Makefile was still buggy, trying to test before the binary had been built. Fixed it and added a diagnostic line to the python script.

4139. By Richard Ferrier on 2014-04-08

Refactored python code in preparation for tests based on the method of manufactured solutions (MMS). Previously the code was littered with nested for-loops to iterate over various simulation options, error norms and convergence rates. This situation onlys worsens with MMS-based tests. The current revision therefore replaces the loops with handling trees that can be composed at run-time to suit the context. The trees are passed command objects encapsulating operations at various levels. These Handler classes and abstract Command class have been placed in test_tools.py. The classes are used in buckley_leverett_test_tools.py (replacing darcy_impes_tools.py) and will be used in manufactured_solution_test_tools.py.

4140. By Richard Ferrier on 2014-04-08

Corrected commented-out line that was breaking a test

4141. By Richard Ferrier on 2014-04-14

Added manufactured_solution_test_tools.py; split up test_tools.py to reflect general-purpose and problem-specific classes

4142. By Richard Ferrier on 2014-04-17

Added the first MMS-based test (the harness will be added shortly). Parameters may be adjusted in due course.

4143. By Richard Ferrier on 2014-04-17

Added code for generating the XML file and cleaning up.

4144. By Richard Ferrier on 2014-04-17

Forgot to add geometry templates to the previous two revisions; these are now in place. Fixed some scripting bugs that were breaking tests.

4145. By Richard Ferrier on 2014-04-19

Geometry templates really checked in this time\!

4146. By Richard Ferrier on 2014-04-23

Minor changes, adding functionality to support factorial experiments where physical variables are controlled.

4147. By Richard Ferrier on 2014-04-23

Continuing the changes of the last revision, generation of symbolic expressions is now done within ManufacturedSolutionTestSuite. Also, WriteXMLSnippet has been replaced by WriteXMLFile to allow multiple solution-dependent XML files to be generated.

4148. By Richard Ferrier on 2014-04-24

Fixed XML inconsistency with schema. Adjusted symbolic expressions to work when generated by a different version of sympy (this generation is not triggered by the test harness, however).

4149. By Richard Ferrier on 2014-04-25

Tweak to the Makefile in darcy_impes_mms_p1_2phase_quadraticrelperm

4150. By Richard Ferrier on 2014-05-01

Added another MMS-based test.

4151. By Richard Ferrier on 2014-05-06

Edited the tests to work around a problem with spatial convergence of the pressure field. The nonconvergence is caused by the pressure fields reconfiguring themselves to be completely different to the forced solution at certain levels of mesh and timestep refinement. This effect has appeared since introducing more timesteps into the tests to make them more rigorous. Saturation fields appear to be unaffected, and they will continue to be used in checking convergence. Note the necessary change of minus sign on diff(phi*self.s[i], t) in darcy_impes_common/manufactured_solution_generation_tools.py. The now-transient tests break otherwise.

4152. By Richard Ferrier on 2014-05-06

Merge with fluidity trunk.

4153. By Richard Ferrier on 2014-05-30

Minor changes.

4154. By Richard Ferrier on 2014-05-30

Merge with fluidity trunk.

Unmerged revisions

4154. By Richard Ferrier on 2014-05-30

Merge with fluidity trunk.

4153. By Richard Ferrier on 2014-05-30

Minor changes.

4152. By Richard Ferrier on 2014-05-06

Merge with fluidity trunk.

4151. By Richard Ferrier on 2014-05-06

Edited the tests to work around a problem with spatial convergence of the pressure field. The nonconvergence is caused by the pressure fields reconfiguring themselves to be completely different to the forced solution at certain levels of mesh and timestep refinement. This effect has appeared since introducing more timesteps into the tests to make them more rigorous. Saturation fields appear to be unaffected, and they will continue to be used in checking convergence. Note the necessary change of minus sign on diff(phi*self.s[i], t) in darcy_impes_common/manufactured_solution_generation_tools.py. The now-transient tests break otherwise.

4150. By Richard Ferrier on 2014-05-01

Added another MMS-based test.

4149. By Richard Ferrier on 2014-04-25

Tweak to the Makefile in darcy_impes_mms_p1_2phase_quadraticrelperm

4148. By Richard Ferrier on 2014-04-24

Fixed XML inconsistency with schema. Adjusted symbolic expressions to work when generated by a different version of sympy (this generation is not triggered by the test harness, however).

4147. By Richard Ferrier on 2014-04-23

Continuing the changes of the last revision, generation of symbolic expressions is now done within ManufacturedSolutionTestSuite. Also, WriteXMLSnippet has been replaced by WriteXMLFile to allow multiple solution-dependent XML files to be generated.

4146. By Richard Ferrier on 2014-04-23

Minor changes, adding functionality to support factorial experiments where physical variables are controlled.

4145. By Richard Ferrier on 2014-04-19

Geometry templates really checked in this time\!

=== modified file 'Makefile.in'
--- Makefile.in	2013-12-13 16:00:27 +0000
+++ Makefile.in	2014-05-30 12:10:35 +0000
@@ -185,6 +185,7 @@
 sw: bin/shallow_water
 be: bin/burgers_equation
 hh: bin/hybridized_helmholtz_solver
+di: bin/darcy_impes
 
 bin/shallow_water: fluidity_library bin/hybridized_helmholtz_solver main/Shallow_Water.F90
 	@cd main; $(MAKE) Shallow_Water.o
@@ -211,6 +212,14 @@
 	@echo "    LD burgers_equation"
 	@$(EVAL) $(FC) -o bin/burgers_equation main/Burgers_Equation.o $(FCFLAGS) $(LIBS)
 
+bin/darcy_impes: fluidity_library main/Darcy_IMPES.F90
+	@cd main; $(MAKE) Darcy_IMPES.o
+	@echo "BUILD darcy_impes"
+	@echo "    MKDIR bin"
+	@mkdir -p bin
+	@echo "    LD darcy_impes"
+	@$(EVAL) $(FC) -o bin/darcy_impes main/Darcy_IMPES.o $(FCFLAGS) $(LIBS)
+
 ifeq (@ENABLE_PSMILE@,yes)
 LIBPSMILE=@FLPSMILE_PATH@/lib/libpsmile_oa4.MPI1.a @FLPSMILE_PATH@/lib/libcommon_oa4.a @FLPSMILE_PATH@/lib/libmpp_io.a -lxml2
 bin/test_coupler: fluidity_library main/test_coupler.F90
@@ -498,12 +507,12 @@
 	@find ./ \( -name make.log \) -exec rm -f {} \; > /dev/null
 	@rm -f Makefile > /dev/null
 
-test: fltools bin/$(FLUIDITY) bin/shallow_water serialtest bin/burgers_equation
+test: fltools bin/$(FLUIDITY) bin/shallow_water bin/darcy_impes serialtest bin/burgers_equation
 
-serialtest: fltools bin/$(FLUIDITY) bin/shallow_water
+serialtest: fltools bin/$(FLUIDITY) bin/shallow_water bin/darcy_impes
 	@cd tests; ../bin/testharness -l short $(EXCLUDE_TAGS) -n $(THREADS)
 
-mediumtest: fltools bin/$(FLUIDITY) manual bin/burgers_equation bin/shallow_water spudtools
+mediumtest: fltools bin/$(FLUIDITY) manual bin/burgers_equation bin/shallow_water bin/darcy_impes spudtools
 	@cd tests; ../bin/testharness -l medium $(EXCLUDE_TAGS) -n $(THREADS)
 
 .PHONY: spudtools
@@ -607,7 +616,7 @@
 	$(MAKE) clean-light
 	@echo "    Congratulations, make makefiles succeeded!"
 
-install: default fltools bin/shallow_water bin/burgers_equation
+install: default fltools bin/shallow_water bin/burgers_equation bin/darcy_impes
 	@mkdir -p $(DESTDIR)$(bindir) $(DESTDIR)$(docdir)/fluidity
 	find bin/ -maxdepth 1 -type f -exec cp '{}' $(DESTDIR)$(bindir) \;
 	rm -f $(DESTDIR)$(bindir)/spud-* $(DESTDIR)$(bindir)/diamond $(DESTDIR)$(bindir)/silenteval $(DESTDIR)$(bindir)/runut

=== added file 'assemble/Darcy_IMPES_Assemble.F90'
--- assemble/Darcy_IMPES_Assemble.F90	1970-01-01 00:00:00 +0000
+++ assemble/Darcy_IMPES_Assemble.F90	2014-05-30 12:10:35 +0000
@@ -0,0 +1,5920 @@
+!    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 darcy_impes_assemble_module
+
+   ! keep in this order, please:
+   use quadrature
+   use elements
+   use sparse_tools
+   use fields
+   use sparse_tools_petsc
+   !
+   use field_derivatives
+   use cv_shape_functions
+   use cv_faces
+   use cvtools
+   use cv_fields
+   use cv_upwind_values
+   use cv_face_values
+   use cv_options
+   use boundary_conditions
+   use fefields, only: compute_cv_mass
+   use petsc_solve_state_module
+   use sparse_matrices_fields
+   use transform_elements, only: transform_cvsurf_to_physical, &
+                                 transform_cvsurf_facet_to_physical
+   use state_module
+   use fldebug
+   use spud
+   use porous_media
+   use parallel_tools
+   use adaptive_timestepping
+   use halos
+   use signal_vars, only : SIG_INT
+   use global_parameters, only : FIELD_NAME_LEN, OPTION_PATH_LEN
+
+   implicit none
+
+   private
+
+   public :: darcy_impes_type, &
+             darcy_impes_cv_options_type, &
+             darcy_impes_copy_to_old, &
+             darcy_impes_copy_to_iterated, &
+             darcy_impes_assemble_and_solve_part_one, &
+             darcy_impes_assemble_and_solve_part_two, &
+             darcy_impes_assemble_and_solve_part_three, &             
+             darcy_impes_calculate_gradient_pressures, &
+             darcy_impes_calculate_non_first_phase_pressures, &
+             darcy_impes_calculate_phase_one_saturation_diagnostic, &
+             darcy_impes_calculate_vel_mob_ff_and_cfl_fields, &
+             darcy_impes_calculate_sum_saturation, &
+             darcy_impes_calculate_relperm_fields, &
+             darcy_impes_calculate_densities, &
+             darcy_impes_calculate_cflnumber_field_based_dt, &
+             darcy_impes_initialise_cached_face_value, &
+             darcy_impes_calculate_inverse_characteristic_length, &
+             darcy_impes_calculate_relperm_den_first_face_values, &
+             darcy_impes_calculate_relperm_isub_face_values, &
+             darcy_impes_get_entire_sfield_boundary_condition, &
+             darcy_impes_get_v_boundary_condition, &
+             darcy_impes_assemble_saturation_rhs_adv, &
+             darcy_impes_assemble_and_solve_phase_pressures, &
+             darcy_impes_calculate_divergence_total_darcy_velocity, &
+             darcy_impes_calculate_inverse_cv_sa
+   
+   ! Parameters defining Darcy IMPES cached options
+   integer, parameter, public :: RELPERM_CORRELATION_POWER                 = 1, &
+                                 RELPERM_CORRELATION_COREY2PHASE           = 2, &
+                                 RELPERM_CORRELATION_COREY2PHASEOPPOSITE   = 3, &
+                                 RELPERM_CORRELATION_MINERAL               = 4, &
+                                 RELPERM_CORRELATION_VANGENUCHTEN          = 5, &
+                                 RELPERM_CORRELATION_JACKSON2PHASE         = 6, &
+                                 RELPERM_CORRELATION_JACKSON2PHASEOPPOSITE = 7
+   
+   ! Parameters defining Darcy IMPES CV face value schemes
+   integer, parameter, public :: DARCY_IMPES_CV_FACEVALUE_NONE                        = 0, &
+                                 DARCY_IMPES_CV_FACEVALUE_FIRSTORDERUPWIND            = 1, &
+                                 DARCY_IMPES_CV_FACEVALUE_FINITEELEMENT               = 2, &
+                                 DARCY_IMPES_CV_FACEVALUE_RELPERMOVERSATUPWIND        = 3, &
+                                 DARCY_IMPES_CV_FACEVALUE_RELPERMOVERSATFINITEELEMENT = 4, &
+                                 DARCY_IMPES_CV_FACEVALUE_RELPERMOVERSATCORRELATION   = 5
+   
+   ! Options associated with the relative permeability correlation
+   type darcy_impes_relperm_corr_options_type
+      ! The correlation type to use
+      integer :: type
+      ! The exponent for the power law correlation for each phase
+      real, dimension(:), pointer :: exponents
+      ! The residual saturation for each phase
+      real, dimension(:), pointer :: residual_saturations
+      ! The cut off value for saturation when calculating relperm
+      real, dimension(:), pointer :: cutoff_saturations
+      ! The scaling_coefficients when calculating relperm
+      real, dimension(:), pointer :: scaling_coefficients      
+   end type darcy_impes_relperm_corr_options_type
+   
+   ! Options associated with the CV discretisation for the darcy impes solver
+   type darcy_impes_cv_options_type
+      ! what CV_FACEVALUE_ to use
+      integer :: facevalue
+      ! the number of face value iterations
+      integer :: number_face_value_iteration
+      ! whether to limit the face value spatially
+      logical :: limit_facevalue
+      ! what CV_LIMITER_ to use
+      integer :: limiter
+      ! the slopes of the limiter being used (if Sweby)
+      real, dimension(2) :: limiter_slopes
+      ! what upwind scheme are we using (if any)
+      integer :: upwind_scheme   
+   end type darcy_impes_cv_options_type
+   
+   ! Options associated with explicit advection subcycling for CV scalar field solver
+   type darcy_impes_subcycle_options_type
+      logical :: have
+      integer :: number
+      logical :: consistent
+   end type darcy_impes_subcycle_options_type
+   
+   ! Options associated with adaptive time stepping
+   type darcy_impes_adaptive_dt_options_type
+      logical :: have
+      real    :: requested_cfl
+      real    :: min_dt
+      real    :: max_dt
+      real    :: increase_tolerance
+      logical :: min_dt_terminate_if_reached
+      logical :: at_first_dt
+   end type darcy_impes_adaptive_dt_options_type
+   
+   ! Options associated with EoS for a phase
+   type darcy_impes_eos_options_type
+      logical :: have_fluids_linear
+      real    :: fluids_linear_reference_density
+   end type darcy_impes_eos_options_type
+   
+   ! Data associated with the cached CV face values
+   type cached_face_value_type
+      real, dimension(:,:,:,:), pointer :: relperm      ! nsub, ngi_vele, vele, nphase
+      real, dimension(:,:,:,:), pointer :: relperm_bdy  ! nsub, ngi_sele, sele, nphase
+      real, dimension(:,:,:),   pointer :: den          !       ngi_vele, vele, nphase
+      real, dimension(:,:,:),   pointer :: den_bdy      !       ngi_sele, sele, nphase
+      real, dimension(:,:),     pointer :: detwei       !       ngi_vele, vele
+      real, dimension(:,:),     pointer :: detwei_bdy   !       ngi_sele, sele
+      real, dimension(:,:,:),   pointer :: normal       ! ndim, ngi_vele, vele
+      real, dimension(:,:,:),   pointer :: normal_bdy   ! ndim, ngi_sele, sele
+      real, dimension(:,:,:,:), pointer :: p_dshape     ! nloc, ngi_vele, ndim, vele
+!!!!! *** THIS IS NOT POSSIBLE YET ***
+!!!!!      real, dimension(:,:,:,:), pointer :: p_dshape_bdy ! nloc, ngi_sele, ndim, sele
+      logical                           :: cached_detwei_normal
+      logical                           :: cached_p_dshape
+   end type cached_face_value_type
+   
+   ! Data associated with a generic prognostic sfield, some of which will point to fields in state
+   type darcy_impes_generic_prog_sfield_type
+      ! The phase this sfield is associated with
+      integer :: phase
+      type(scalar_field), pointer :: sfield
+      type(scalar_field), pointer :: old_sfield
+      type(scalar_field), pointer :: sfield_abs
+      type(scalar_field), pointer :: sfield_src
+      type(tensor_field), pointer :: sfield_diff
+      logical :: have_diff
+      logical :: have_abs
+      logical :: have_src
+      logical :: have_adv
+      type(darcy_impes_cv_options_type) :: sfield_cv_options
+   end type darcy_impes_generic_prog_sfield_type
+   
+   type darcy_impes_type
+      ! *** Pointers to fields from state that have array length of number of phases ***
+      type(vector_field_pointer), dimension(:), pointer :: darcy_velocity
+      type(scalar_field_pointer), dimension(:), pointer :: mobility      
+      type(scalar_field_pointer), dimension(:), pointer :: fractional_flow
+      type(scalar_field_pointer), dimension(:), pointer :: saturation
+      type(scalar_field_pointer), dimension(:), pointer :: old_saturation
+      type(scalar_field_pointer), dimension(:), pointer :: saturation_source
+      type(scalar_field_pointer), dimension(:), pointer :: relative_permeability
+      type(scalar_field_pointer), dimension(:), pointer :: old_relative_permeability
+      type(scalar_field_pointer), dimension(:), pointer :: viscosity
+      type(scalar_field_pointer), dimension(:), pointer :: cfl
+      type(scalar_field_pointer), dimension(:), pointer :: pressure
+      type(scalar_field_pointer), dimension(:), pointer :: capilliary_pressure
+      type(scalar_field_pointer), dimension(:), pointer :: density
+      type(scalar_field_pointer), dimension(:), pointer :: old_density
+      ! *** Pointers to fields from state that are NOT phase dependent ***
+      type(mesh_type),    pointer :: pressure_mesh
+      type(mesh_type),    pointer :: elementwise_mesh
+      type(scalar_field), pointer :: average_pressure
+      type(scalar_field), pointer :: porosity
+      type(scalar_field), pointer :: absolute_permeability
+      type(vector_field), pointer :: positions
+      type(vector_field), pointer :: total_darcy_velocity
+      type(scalar_field), pointer :: total_mobility
+      type(scalar_field), pointer :: sum_saturation
+      type(scalar_field), pointer :: div_total_darcy_velocity
+      type(vector_field), pointer :: gravity_direction
+      type(vector_field), pointer :: bulk_darcy_velocity
+      ! ***** Pointers to DUAL fields from state that have array length of number of phases *****
+      type(scalar_field_pointer), dimension(:), pointer :: pressure_other_porous_media
+      type(scalar_field_pointer), dimension(:), pointer :: transmissibility_lambda_dual
+      ! *** Pointer to the pressure mesh - pressure mesh sparsity, used for pressure matrix and finding CV upwind values ***
+      type(csr_sparsity), pointer :: sparsity_pmesh_pmesh
+      ! *** Data associated with generic prognostic scalar fields, some of which points to fields in state ***
+      type(darcy_impes_generic_prog_sfield_type), dimension(:), pointer :: generic_prog_sfield
+      ! *** Data associated with GradientPressure fields ***
+      type(element_type)                                :: gradient_pressure_shape
+      type(mesh_type)                                   :: gradient_pressure_mesh
+      type(vector_field_pointer), dimension(:), pointer :: gradient_pressure
+      type(vector_field_pointer), dimension(:), pointer :: iterated_gradient_pressure
+      ! *** Data associated with gravity ***
+      logical            :: have_gravity
+      real               :: gravity_magnitude
+      type(vector_field) :: gravity
+      ! *** Fields allocated here used in assemble algorithm ***
+      type(vector_field)     :: positions_pressure_mesh
+      type(csr_matrix)       :: matrix
+      type(petsc_csr_matrix) :: dual_block_pressure_matrix
+      type(csr_matrix)       :: pressure_matrix
+      type(scalar_field), pointer :: pressure_rhs
+      type(scalar_field)     :: lhs
+      type(scalar_field)     :: rhs
+      type(scalar_field)     :: rhs_full
+      type(scalar_field)     :: rhs_high_resolution
+      type(scalar_field)     :: inverse_cv_mass_pressure_mesh
+      type(scalar_field)     :: cv_mass_pressure_mesh_with_source
+      type(scalar_field)     :: cv_mass_pressure_mesh_with_porosity   
+      type(scalar_field)     :: cv_mass_pressure_mesh_with_old_porosity 
+      type(scalar_field)     :: inverse_cv_sa_pressure_mesh
+      type(scalar_field)     :: modified_relative_permeability
+      type(scalar_field), pointer :: constant_zero_sfield_pmesh
+      type(scalar_field), pointer :: constant_zero_sfield_elementwisemesh
+      type(scalar_field_pointer), dimension(:), pointer :: old_saturation_subcycle
+      ! ***** coupling field associated with DUAL model *****
+      type(scalar_field) :: cv_mass_pressure_mesh_with_lambda_dual
+      ! ***** rhs src field associated with DUAL model *****
+      type(scalar_field) :: rhs_dual
+      ! *** Data associated with v, pressure and sfield BC allocated here ***
+      type(mesh_type)                          :: bc_surface_mesh
+      type(scalar_field)                       :: v_bc_value
+      integer,           dimension(:), pointer :: v_bc_flag
+      type(scalar_field)                       :: sfield_bc_value
+      integer,           dimension(:), pointer :: sfield_bc_flag
+      type(scalar_field)                       :: pressure_bc_value
+      integer,           dimension(:), pointer :: pressure_bc_flag
+      type(scalar_field)                       :: inverse_characteristic_length
+      real                                     :: weak_pressure_bc_coeff
+      ! *** The number of phase and the CV surface quadrature degree to use ***
+      integer :: number_phase, cv_surface_quaddegree   
+      ! *** Data specifically associated with the CV discretisation allocated here ***
+      type(darcy_impes_cv_options_type)       :: saturation_cv_options
+      type(darcy_impes_cv_options_type)       :: relperm_cv_options
+      type(darcy_impes_cv_options_type)       :: density_cv_options
+      type(cv_faces_type)                     :: cvfaces
+      type(element_type)                      :: x_cvshape_full
+      type(element_type)                      :: p_cvshape_full
+      type(element_type)                      :: gradp_cvshape_full
+      type(element_type)                      :: x_cvshape
+      type(element_type)                      :: p_cvshape
+      type(element_type)                      :: gradp_cvshape
+      type(element_type)                      :: x_cvbdyshape
+      type(element_type)                      :: p_cvbdyshape      
+      type(element_type)                      :: gradp_cvbdyshape
+      type(darcy_impes_subcycle_options_type) :: subcy_opt_sat
+      type(csr_matrix)                        :: relperm_upwind
+      type(csr_matrix)                        :: sfield_upwind
+      ! *** The minimum value of saturation to use in the denominator of modrelperm face value ***
+      real :: minimum_denominator_saturation_value
+      ! *** A flag for whether saturation face values need to be determined
+      logical :: determine_saturation_face_values
+      ! *** Data associate with the relperm correlation options ***
+      type(darcy_impes_relperm_corr_options_type) :: relperm_corr_options
+      ! *** Flag for each phase for whether there is capilliary pressure ***
+      logical, dimension(:), pointer :: have_capilliary_pressure
+      ! *** Flag for each phase for whether there is a saturation source ***
+      logical, dimension(:), pointer :: have_saturation_source
+      ! *** Flag for whether the first phase saturation is diagnostic, else it is prognostic ***
+      logical :: phase_one_saturation_diagnostic      
+      ! *** Flag for whether the first phase pressure is prognostic, else it is prescribed *** 
+      logical :: first_phase_pressure_prognostic
+      ! *** The cached face values for all phases ***
+      type(cached_face_value_type) :: cached_face_value
+      ! *** The advection subcycle timestep size ***
+      real    :: dt_subcycle
+      ! *** Time time step size, also stored here for convenience ***
+      real :: dt
+      ! *** The current time, also stored here for convenience ***
+      real :: current_time 
+      ! *** Non linear iteration, also stored here for convencience ***
+      integer :: nonlinear_iter
+      ! *** Max Non linear iteration for this timestep, also stored here for convencience ***
+      integer :: max_nonlinear_iter_this_timestep
+      ! *** The number of volume, surface elements and nodes, stored here for convenience ***
+      integer :: number_vele
+      integer :: number_sele
+      integer :: number_pmesh_node
+      ! *** Geometric dimension, also stored here for convenience ***
+      integer :: ndim
+      ! *** Options data associated with adaptive time stepping stored here ***
+      type(darcy_impes_adaptive_dt_options_type) :: adaptive_dt_options
+      ! *** Options data associated with each phase EoS stored here *** 
+      type(darcy_impes_eos_options_type), dimension(:), pointer :: eos_options
+      ! *** Pointer to main state array, for convenience ***
+      type(state_type), dimension(:), pointer :: state
+   end type darcy_impes_type
+   
+   contains
+
+! ----------------------------------------------------------------------------
+  
+   subroutine darcy_impes_copy_to_old(di)
+            
+      !!< Copy to old darcy impes data not associated with di%state
+
+      type(darcy_impes_type), intent(inout) :: di
+            
+      call set(di%cv_mass_pressure_mesh_with_old_porosity, di%cv_mass_pressure_mesh_with_porosity)
+      
+   end subroutine darcy_impes_copy_to_old
+
+! ----------------------------------------------------------------------------
+  
+   subroutine darcy_impes_copy_to_iterated(di)
+            
+      !!< Copy to iterated darcy impes data not associated with di%state
+
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variables
+      integer :: p
+      
+      do p = 1, di%number_phase
+      
+         call set(di%iterated_gradient_pressure(p)%ptr, di%gradient_pressure(p)%ptr)
+      
+      end do
+      
+   end subroutine darcy_impes_copy_to_iterated
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_assemble_and_solve_part_one(di, have_dual)
+      
+      !!< Assemble and solve the Darcy equations using an CIMPESS algorithm
+      !!< which is a modification of the IMPES to include consistent subcycling.
+      !!< Part one: form cv mass and initial saturation solve for subcycling
+      
+      type(darcy_impes_type), intent(inout) :: di
+      logical ,               intent(in)    :: have_dual      
+      
+      ! local variables
+      logical :: form_new_subcycle_relperm_face_values
+      
+      ewrite(1,*) 'Start Darcy IMPES assemble and solve part one'
+                  
+      ! Calculate the latest CV mass on the pressure mesh with porosity included
+      call compute_cv_mass(di%positions, di%cv_mass_pressure_mesh_with_porosity, di%porosity)      
+      
+      ! If it is the first iteration and there is subcycling (>1)
+      ! and a consistent global continuity then solve the saturations 
+      if ((di%nonlinear_iter == 1) .and. &
+          (di%subcy_opt_sat%number > 1) .and. &
+          (di%subcy_opt_sat%consistent)) then
+         
+         form_new_subcycle_relperm_face_values  = .true.
+         
+         call darcy_impes_assemble_and_solve_phase_saturations(di, &
+                                                               form_new_subcycle_relperm_face_values, &
+                                                               have_dual)
+         
+      end if
+                        
+      ewrite(1,*) 'Finished Darcy IMPES assemble and solve part one'
+           
+   end subroutine darcy_impes_assemble_and_solve_part_one
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_assemble_and_solve_part_two(di, di_dual, have_dual, solve_dual_pressure)
+      
+      !!< Assemble and solve the Darcy equations using an CIMPESS algorithm
+      !!< which is a modification of the IMPES to include consistent subcycling.
+      !!< Part two: solve for the pressure
+      
+      type(darcy_impes_type), intent(inout) :: di
+      type(darcy_impes_type), intent(inout) :: di_dual
+      logical ,               intent(in)    :: have_dual      
+      logical ,               intent(in)    :: solve_dual_pressure
+            
+      ewrite(1,*) 'Start Darcy IMPES assemble and solve part two'
+                  
+      ! Assemble and solve the phase pressures
+      call darcy_impes_assemble_and_solve_phase_pressures(di, di_dual, have_dual, solve_dual_pressure)
+                  
+      ewrite(1,*) 'Finished Darcy IMPES assemble and solve part two'
+           
+   end subroutine darcy_impes_assemble_and_solve_part_two
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_assemble_and_solve_part_three(di, have_dual)
+      
+      !!< Assemble and solve the Darcy equations using an CIMPESS algorithm
+      !!< which is a modification of the IMPES to include consistent subcycling.
+      !!< Part three: find pressure gradients, calc div.u, solve for saturation
+      !!<             and generic scalar fields, calc sum S, calc density, calc relperm.
+      
+      type(darcy_impes_type), intent(inout) :: di
+      logical ,               intent(in)    :: have_dual      
+      
+      ! local variables
+      logical :: form_new_subcycle_relperm_face_values
+      
+      ewrite(1,*) 'Start Darcy IMPES assemble and solve part three'
+                  
+      ! Calculate the gradient pressures 
+      call darcy_impes_calculate_gradient_pressures(di)
+
+      ! Calculate the divergence total darcy velocity
+      call darcy_impes_calculate_divergence_total_darcy_velocity(di)
+      
+      ! Assemble and solve the phase saturations
+            
+      if ((di%nonlinear_iter == di%max_nonlinear_iter_this_timestep) .and. &
+          (di%subcy_opt_sat%consistent)) then
+
+         form_new_subcycle_relperm_face_values = .false.
+
+      else 
+
+         form_new_subcycle_relperm_face_values = .true.
+
+      end if 
+            
+      call darcy_impes_assemble_and_solve_phase_saturations(di, &
+                                                            form_new_subcycle_relperm_face_values, &
+                                                            have_dual)
+      
+      call darcy_impes_assemble_and_solve_generic_prog_sfields(di)
+      
+      ! Calculate the sum of the saturations
+      call darcy_impes_calculate_sum_saturation(di)
+      
+      ! Calculate the density field of each phase
+      call darcy_impes_calculate_densities(di)
+
+      ! Calculate the relative permeabilities of each phase
+      call darcy_impes_calculate_relperm_fields(di)
+            
+      ewrite(1,*) 'Finished Darcy IMPES assemble and solve part three'
+           
+   end subroutine darcy_impes_assemble_and_solve_part_three
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_assemble_and_solve_phase_pressures(di, di_dual, have_dual, solve_dual_pressure)
+      
+      !!< Assemble and solve the phase pressures. The first phase pressure 
+      !!< is solved via a matrix equation (if prognostic) and the other
+      !!< phases pressures are calculated from the first phase pressure
+      !!< and their own capilliary pressures
+      
+      type(darcy_impes_type), intent(inout) :: di
+      type(darcy_impes_type), intent(inout) :: di_dual
+      logical ,               intent(in)    :: have_dual  
+      logical ,               intent(in)    :: solve_dual_pressure
+      
+      ! Assemble and solve the first phase pressure if it is prognostic
+      if (di%first_phase_pressure_prognostic) call darcy_impes_assemble_and_solve_first_phase_pressure(di, di_dual, have_dual, solve_dual_pressure)
+      
+      ! Calculate the non first phase pressure's
+      call darcy_impes_calculate_non_first_phase_pressures(di) 
+      if (have_dual)  call darcy_impes_calculate_non_first_phase_pressures(di_dual)     
+      
+      ! Calculate the average pressure
+      call darcy_impes_calculate_average_pressure(di)
+      if (have_dual) call darcy_impes_calculate_average_pressure(di_dual)
+      
+   end subroutine darcy_impes_assemble_and_solve_phase_pressures
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_assemble_and_solve_first_phase_pressure(di, di_dual, have_dual, solve_dual_pressure)
+      
+      !!< Assemble and solve the first phase pressure. A source is included 
+      !!< due to the capilliary pressures of non first phases as well as gravity 
+      !!< and the rate of change of porosity and the individual phase sources.
+      !!< Face values for relative permeability and density have already been evaluated.
+      !!< If have_dual and solve_dual_pressure are both true then a block matrix using 
+      !!< di%dual_block_pressure_matrix is solved for. If have_dual is true but 
+      !!< solve_dual_pressure is false then a single matrix using di%pressure_matrix 
+      !!< is solved for the prime pressure phase 1 then the dual pressure phase 1 
+      !!< is set to this. If have_dual is false then di%pressure_matrix is used.
+      !!< (This is because an error is occuring for petsc_csr_matrix in 2d)
+      
+      type(darcy_impes_type), intent(inout) :: di
+      type(darcy_impes_type), intent(inout) :: di_dual
+      logical ,               intent(in)    :: have_dual  
+      logical ,               intent(in)    :: solve_dual_pressure
+      
+      ! local variables
+      type(scalar_field_pointer), dimension(:), pointer :: all_pressure
+      type(scalar_field_pointer), dimension(:), pointer :: all_rhs
+      integer :: p, vele, sele, iloc, oloc, jloc, face, gi, ggi
+      real    :: face_value, grad_cap_p_dot_n, g_dot_n
+      real,    dimension(1)              :: absperm_ele, visc_ele
+      real,    dimension(:,:),   pointer :: grav_ele
+      real,    dimension(:,:),   pointer :: grad_cap_pressure_face_quad
+      real,    dimension(:,:),   pointer :: x_ele
+      real,    dimension(:,:,:), pointer :: p_dshape
+      real,    dimension(:,:),   pointer :: normal
+      real,    dimension(:),     pointer :: detwei
+      real,    dimension(:),     pointer :: normgi
+      logical, dimension(:),     pointer :: notvisited
+      real,    dimension(:,:),   pointer :: p_mat_local
+      real,    dimension(:),     pointer :: p_rhs_local
+      real,    dimension(:,:),   pointer :: x_face_quad
+      integer, dimension(:),     pointer :: x_pmesh_nodes
+      integer, dimension(:),     pointer :: p_nodes      
+      real,    dimension(1)              :: absperm_ele_bdy, visc_ele_bdy
+      real,    dimension(:,:),   pointer :: grav_ele_bdy
+!!!!! *** THIS IS NOT POSSIBLE YET ***
+!!!!!      real,    dimension(:,:),   pointer :: grad_cap_pressure_face_quad_bdy
+!!!!!      real,    dimension(:,:,:), pointer :: p_dshape_bdy
+      real,    dimension(:),     pointer :: bc_sele_val
+      real,    dimension(:),     pointer :: inv_char_len_ele_bdy
+      real,    dimension(:,:),   pointer :: normal_bdy
+      real,    dimension(:),     pointer :: detwei_bdy
+      real,    dimension(:,:),   pointer :: x_ele_bdy
+      real,    dimension(:),     pointer :: p_rhs_local_bdy
+      real,    dimension(:,:),   pointer :: p_matrix_local_bdy
+      integer, dimension(:),     pointer :: p_nodes_bdy
+      integer, parameter :: PRESSURE_BC_TYPE_WEAKDIRICHLET   = 1, PRESSURE_BC_TYPE_DIRICHLET = 2 
+      integer, parameter :: V_BC_TYPE_PRESCRIBED_NORMAL_FLOW = 1, V_BC_TYPE_NO_NORMAL_FLOW   = 2
+            
+      ! allocate arrays used in assemble process - many assume
+      ! that all elements are the same type.
+      allocate(x_ele(di%ndim, ele_loc(di%positions,1)))      
+      if (.not. di%cached_face_value%cached_p_dshape) then
+         allocate(p_dshape(ele_loc(di%pressure_mesh,1), di%p_cvshape%ngi, mesh_dim(di%pressure_mesh)))
+      end if
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(normal(di%ndim,di%x_cvshape%ngi))
+         allocate(detwei(di%x_cvshape%ngi))      
+      end if
+      allocate(normgi(di%ndim))
+      allocate(notvisited(di%x_cvshape%ngi))
+      allocate(p_mat_local(ele_loc(di%pressure_mesh,1), ele_loc(di%pressure_mesh,1)))
+      allocate(p_rhs_local(ele_loc(di%pressure_mesh,1)))
+      allocate(x_face_quad(di%ndim, di%x_cvshape%ngi))
+      allocate(grav_ele(di%ndim,1))
+      allocate(grad_cap_pressure_face_quad(di%ndim, di%p_cvshape%ngi))
+
+      allocate(grav_ele_bdy(di%ndim,1))
+      allocate(bc_sele_val(face_loc(di%pressure_mesh,1)))
+      allocate(inv_char_len_ele_bdy(face_loc(di%pressure_mesh,1)))
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(detwei_bdy(di%x_cvbdyshape%ngi))
+         allocate(normal_bdy(di%ndim, di%x_cvbdyshape%ngi))
+      end if
+      allocate(p_rhs_local_bdy(face_loc(di%pressure_mesh,1)))
+      allocate(p_matrix_local_bdy(face_loc(di%pressure_mesh,1),face_loc(di%pressure_mesh,1)))
+      allocate(x_ele_bdy(di%ndim, face_loc(di%positions,1)))      
+      allocate(p_nodes_bdy(face_loc(di%pressure_mesh,1)))
+!!!!! *** THIS IS NOT POSSIBLE YET ***
+!!!!!      allocate(grad_cap_pressure_face_quad_bdy(di%ndim, di%p_cvbdyshape_full%ngi))
+!!!!!      if (.not. di%cached_face_value%cached_p_dshape) then
+!!!!!         allocate(p_dshape_bdy(face_loc(di%pressure_mesh,1), di%p_cvbdyshape_full%ngi, mesh_dim(di%pressure_mesh)))
+!!!!!      end if
+      
+      ewrite(1,*) 'Solve first phase Pressure'
+      
+      ! Initialise the matrix and rhs and solution pointer
+      call zero(di%pressure_rhs)      
+      
+      if (have_dual .and. solve_dual_pressure) then
+         call zero(di%dual_block_pressure_matrix)
+         call zero(di_dual%pressure_rhs)
+         
+         allocate(all_pressure(2))
+         allocate(all_rhs(2))
+         all_pressure(1)%ptr => di%pressure(1)%ptr
+         all_pressure(2)%ptr => di_dual%pressure(1)%ptr
+         all_rhs(1)%ptr => di%pressure_rhs
+         all_rhs(2)%ptr => di_dual%pressure_rhs         
+      else
+         call zero(di%pressure_matrix)
+         allocate(all_pressure(1))
+         allocate(all_rhs(1)) 
+         all_pressure(1)%ptr => di%pressure(1)%ptr
+         all_rhs(1)%ptr => di%pressure_rhs                 
+      end if
+
+      ewrite(1,*) 'Add phase sources to global continuity'
+      
+      src_phase_loop: do p = 1, di%number_phase
+         
+         ! Should this take account of subcycling?!?!
+         call compute_cv_mass(di%positions, di%cv_mass_pressure_mesh_with_source, di%saturation_source(p)%ptr)
+         
+         ! Should this include porosity ?!?!
+         call addto(di%pressure_rhs, di%cv_mass_pressure_mesh_with_source)
+      
+      end do src_phase_loop
+      
+      if (have_dual .and. solve_dual_pressure) then
+      
+         dual_src_phase_loop: do p = 1, di_dual%number_phase
+
+            ! Should this take account of subcycling?!?!
+            call compute_cv_mass(di_dual%positions, di_dual%cv_mass_pressure_mesh_with_source, di_dual%saturation_source(p)%ptr)
+
+            ! Should this include porosity ?!?!
+            call addto(di_dual%pressure_rhs, di_dual%cv_mass_pressure_mesh_with_source)
+
+         end do dual_src_phase_loop
+      
+      end if
+
+      ! Add the transmissibility lambda dual to each block diagonal      
+      if (have_dual .and. solve_dual_pressure) then         
+         transmiss_phase_loop: do p = 1, di_dual%number_phase
+            call compute_cv_mass(di_dual%positions, di_dual%cv_mass_pressure_mesh_with_lambda_dual, di_dual%transmissibility_lambda_dual(p)%ptr)
+
+            call addto_diag(di%dual_block_pressure_matrix, 1, 1, di_dual%cv_mass_pressure_mesh_with_lambda_dual, scale =  1.0)
+            call addto_diag(di%dual_block_pressure_matrix, 1, 2, di_dual%cv_mass_pressure_mesh_with_lambda_dual, scale = -1.0)
+            call addto_diag(di%dual_block_pressure_matrix, 2, 1, di_dual%cv_mass_pressure_mesh_with_lambda_dual, scale = -1.0)
+            call addto_diag(di%dual_block_pressure_matrix, 2, 2, di_dual%cv_mass_pressure_mesh_with_lambda_dual, scale =  1.0)
+            
+            if (p > 1) then
+               ! Add lambda*(PC_dual - PC_prime) to prime rhs
+               call set(di%rhs, di_dual%capilliary_pressure(p)%ptr)
+               call addto(di%rhs, di%capilliary_pressure(p)%ptr, scale = -1.0)
+               call scale(di%rhs, di_dual%cv_mass_pressure_mesh_with_lambda_dual)               
+               call addto(di%pressure_rhs, di%rhs)
+
+               ! Add lambda*(PC_prime- PC_dual) to dual rhs
+               call set(di%rhs, di%capilliary_pressure(p)%ptr)
+               call addto(di%rhs, di_dual%capilliary_pressure(p)%ptr, scale = -1.0)
+               call scale(di%rhs, di_dual%cv_mass_pressure_mesh_with_lambda_dual)
+               
+               call addto(di_dual%pressure_rhs, di%rhs)               
+            end if
+         end do transmiss_phase_loop
+      end if
+                  
+      ewrite(1,*) 'Add rate of change of porosity to global continuity equation'
+      
+      ! Add rate of change of porosity to rhs    
+      call addto(di%pressure_rhs, di%cv_mass_pressure_mesh_with_old_porosity, scale = 1.0/di%dt)
+                  
+      call addto(di%pressure_rhs, di%cv_mass_pressure_mesh_with_porosity, scale = -1.0/di%dt)
+      
+      if (have_dual .and. solve_dual_pressure) then
+      
+         call addto(di_dual%pressure_rhs, di_dual%cv_mass_pressure_mesh_with_old_porosity, scale = 1.0/di%dt)
+                  
+         call addto(di_dual%pressure_rhs, di_dual%cv_mass_pressure_mesh_with_porosity, scale = -1.0/di%dt)
+      
+      end if
+      
+      ! Assemble a contribution from each phase to form a global continuity equation to solve for first phase pressure
+      phase_loop: do p = 1, di%number_phase
+         
+         ewrite(1,*) 'Assemble volume contribution to global continuity from phase ',p
+            
+         ! Loop volume elements assembling local contributions     
+         vol_element_loop: do vele = 1, di%number_vele
+                       
+            call pressure_volume_element_local_assemble(di) 
+            
+            ! Add volume element contribution to global pressure matrix and rhs            
+            if (have_dual .and. solve_dual_pressure) then 
+               call addto(di%dual_block_pressure_matrix, 1, 1, p_nodes, p_nodes, p_mat_local)
+            
+               call addto(di%pressure_rhs, p_nodes, p_rhs_local)
+ 
+               call pressure_volume_element_local_assemble(di_dual) 
+               
+               call addto(di%dual_block_pressure_matrix, 2, 2, p_nodes, p_nodes, p_mat_local)
+               
+               call addto(di_dual%pressure_rhs, p_nodes, p_rhs_local)
+            
+            else 
+
+              call addto(di%pressure_matrix, p_nodes, p_nodes, p_mat_local)
+            
+              call addto(di%pressure_rhs, p_nodes, p_rhs_local)
+            
+            end if
+
+         end do vol_element_loop
+         
+      end do phase_loop
+                  
+      ! Add normal v BC integrals for each phase, else include if required 
+      ! a weak pressure BC (which if not given is assumed zero).
+
+      ! Get the pressure BC - if no v and weak pressure then extra integrals are added
+      ! - also get strong pressure BC as this implies no need to add any surface integrals
+      call darcy_impes_get_entire_sfield_boundary_condition(di%pressure(1)%ptr, &
+                                                            (/"weakdirichlet", &
+                                                              "dirichlet    "/), &
+                                                            di%bc_surface_mesh, &
+                                                            di%pressure_bc_value, &
+                                                            di%pressure_bc_flag)
+      
+      if (have_dual .and. solve_dual_pressure) then
+         call darcy_impes_get_entire_sfield_boundary_condition(di_dual%pressure(1)%ptr, &
+                                                               (/"weakdirichlet", &
+                                                                 "dirichlet    "/), &
+                                                               di_dual%bc_surface_mesh, &
+                                                               di_dual%pressure_bc_value, &
+                                                               di_dual%pressure_bc_flag)      
+      end if 
+      
+      phase_loop_bc: do p = 1, di%number_phase
+
+         ewrite(1,*) 'Assemble boundary contribution to global continuity from phase ',p
+                           
+         ! Get this phase v BC info - only for no_normal_flow and prescribed_normal_flow which is special as it is a scalar
+         call darcy_impes_get_v_boundary_condition(di%darcy_velocity(p)%ptr, &
+                                                   (/"prescribed_normal_flow", &
+                                                     "no_normal_flow        "/), &
+                                                   di%bc_surface_mesh, &
+                                                   di%v_bc_value, &
+                                                   di%v_bc_flag)
+         
+         if (have_dual .and. solve_dual_pressure) then         
+            call darcy_impes_get_v_boundary_condition(di_dual%darcy_velocity(p)%ptr, &
+                                                      (/"prescribed_normal_flow", &
+                                                        "no_normal_flow        "/), &
+                                                      di_dual%bc_surface_mesh, &
+                                                      di_dual%v_bc_value, &
+                                                      di_dual%v_bc_flag)
+         end if
+
+         sele_loop: do sele = 1, di%number_sele
+            
+            call pressure_surface_element_local_assemble(di)
+            
+            if (have_dual .and. solve_dual_pressure) then            
+
+               call addto(di%dual_block_pressure_matrix, 1, 1, p_nodes_bdy, p_nodes_bdy, p_matrix_local_bdy)
+
+               call addto(di%pressure_rhs, p_nodes_bdy, p_rhs_local_bdy)
+
+               call pressure_surface_element_local_assemble(di_dual) 
+                       
+               call addto(di%dual_block_pressure_matrix, 2, 2, p_nodes_bdy, p_nodes_bdy, p_matrix_local_bdy)
+
+               call addto(di_dual%pressure_rhs, p_nodes_bdy, p_rhs_local_bdy)            
+            
+            else
+
+               call addto(di%pressure_matrix, p_nodes_bdy, p_nodes_bdy, p_matrix_local_bdy)
+
+               call addto(di%pressure_rhs, p_nodes_bdy, p_rhs_local_bdy)
+            
+            end if
+
+         end do sele_loop      
+
+      end do phase_loop_bc
+            
+      ! Apply any strong dirichlet BC's - that can only be applied to the first phase pressure
+      ! of both the prime and dual porous media
+      if (have_dual .and. solve_dual_pressure) then
+         call apply_dirichlet_conditions(di%dual_block_pressure_matrix, all_rhs, all_pressure)
+      else
+         call apply_dirichlet_conditions(di%pressure_matrix, all_rhs(1)%ptr, all_pressure(1)%ptr)      
+      end if
+      
+      ! Solve the pressure(s)
+      if (have_dual .and. solve_dual_pressure) then
+         call petsc_solve(all_pressure, di%dual_block_pressure_matrix, all_rhs, di%state(1))
+      else
+         call petsc_solve(all_pressure(1)%ptr, di%pressure_matrix, all_rhs(1)%ptr, di%state(1))      
+      end if
+      
+      ! Set the strong BC nodes to the values to be consistent
+      call set_dirichlet_consistent(di%pressure(1)%ptr) 
+      
+      if (have_dual) then
+         
+         if (solve_dual_pressure) then
+         
+            call set_dirichlet_consistent(di_dual%pressure(1)%ptr) 
+         
+         else
+            
+            call set(di_dual%pressure(1)%ptr, di%pressure(1)%ptr)
+         
+         end if
+         
+      end if
+          
+      ! deallocate local variables as required
+      deallocate(x_ele)
+      if (di%cached_face_value%cached_p_dshape) then
+         nullify(p_dshape)      
+      else
+         deallocate(p_dshape)
+      end if
+      if (di%cached_face_value%cached_detwei_normal) then
+         nullify(detwei)
+         nullify(normal)      
+      else
+         deallocate(detwei)
+         deallocate(normal)
+      end if
+      deallocate(normgi)
+      deallocate(notvisited)
+      deallocate(p_mat_local)
+      deallocate(p_rhs_local)
+      deallocate(x_face_quad)
+      deallocate(grav_ele)
+      deallocate(grad_cap_pressure_face_quad)
+      
+      deallocate(grav_ele_bdy)
+      deallocate(bc_sele_val)
+      deallocate(inv_char_len_ele_bdy)
+      if (di%cached_face_value%cached_detwei_normal) then
+         nullify(detwei_bdy)
+         nullify(normal_bdy)      
+      else
+         deallocate(detwei_bdy)
+         deallocate(normal_bdy)
+      end if
+      deallocate(p_rhs_local_bdy)
+      deallocate(p_matrix_local_bdy)
+      deallocate(x_ele_bdy)      
+      deallocate(p_nodes_bdy)
+!!!!! *** THIS IS NOT POSSIBLE YET ***
+!!!!!      deallocate(grad_cap_pressure_face_quad_bdy)
+
+      deallocate(all_pressure)
+      deallocate(all_rhs)
+      
+      ewrite(1,*) 'Finished solve first phase Pressure'
+      
+      ewrite_minmax(di%pressure(1)%ptr)
+      
+      if (have_dual) then
+         ewrite_minmax(di_dual%pressure(1)%ptr)
+      end if
+      
+      contains
+      
+      ! --------------------------------------------------------------------------------------
+      
+         subroutine pressure_volume_element_local_assemble(di) 
+
+            !!< Assemble the local volume element pressure matrix and rhs
+         
+            type(darcy_impes_type), intent(inout) :: di
+                        
+            ! The node indices of the pressure mesh
+            p_nodes => ele_nodes(di%pressure_mesh, vele)
+            
+            ! The gravity values for this element for each direction
+            grav_ele = ele_val(di%gravity, vele) 
+
+            ! get the viscosity value for this element for this phase
+            visc_ele = ele_val(di%viscosity(p)%ptr, vele)
+
+            ! get the absolute permeability value for this element
+            absperm_ele = ele_val(di%absolute_permeability, vele)         
+            
+            ! obtain the transformed determinant*weight and normals
+            if (di%cached_face_value%cached_detwei_normal) then
+               
+               detwei => di%cached_face_value%detwei(:,vele)
+               
+               normal => di%cached_face_value%normal(:,:,vele)
+               
+            else
+
+               ! get the coordinate values for this element for each positions local node
+               x_ele = ele_val(di%positions, vele)         
+
+               ! get the coordinate values for this element for each quadrature point
+               x_face_quad = ele_val_at_quad(di%positions, vele, di%x_cvshape)
+
+               ! The node indices of the positions projected to the pressure mesh
+               x_pmesh_nodes => ele_nodes(di%positions_pressure_mesh, vele)
+            
+               call transform_cvsurf_to_physical(x_ele, di%x_cvshape, detwei, normal, di%cvfaces)
+            
+            end if
+            
+            ! obtain the derivative of the pressure mesh shape function at the CV face quadrature points
+            if (di%cached_face_value%cached_p_dshape) then
+               
+               p_dshape => di%cached_face_value%p_dshape(:,:,:,vele)
+               
+            else
+            
+               call transform_to_physical(di%positions, vele, x_shape = di%x_cvshape_full, &
+                                          shape = di%p_cvshape_full, dshape = p_dshape)
+            
+            end if
+            
+            ! get the gradient capilliary pressure at the cv surface quadrature points for each direction 
+            if (p > 1) then
+               grad_cap_pressure_face_quad = &
+              &darcy_impes_ele_grad_at_quad_scalar(di%capilliary_pressure(p)%ptr, vele, dn = p_dshape)
+            end if
+            
+            ! Initialise array for the quadrature points of this 
+            ! element for whether it has already been visited
+            notvisited = .true.
+
+            ! Initialise the local p matrix and rhs to assemble for this element
+            p_mat_local = 0.0
+            p_rhs_local = 0.0
+            
+            ! loop over local nodes within this element
+            nodal_loop_i: do iloc = 1, di%pressure_mesh%shape%loc
+
+               ! loop over CV faces internal to this element
+               face_loop: do face = 1, di%cvfaces%faces
+
+                  ! is this a face neighbouring iloc?
+                  is_neigh: if(di%cvfaces%neiloc(iloc, face) /= 0) then
+
+                     ! find the opposing local node across the CV face
+                     oloc = di%cvfaces%neiloc(iloc, face)
+
+                     ! loop over gauss points on face
+                     quadrature_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                        ! global gauss pt index for this element
+                        ggi = (face-1)*di%cvfaces%shape%ngi + gi
+
+                        ! check if this quadrature point has already been visited
+                        check_visited: if(notvisited(ggi)) then
+
+                           notvisited(ggi) = .false.
+
+                           if (di%cached_face_value%cached_detwei_normal) then
+
+                              ! cached normal has correct orientation already
+                              normgi = normal(:,ggi)
+
+                           else
+
+                              ! correct the orientation of the normal so it points away from iloc
+                              normgi = orientate_cvsurf_normgi(node_val(di%positions_pressure_mesh, x_pmesh_nodes(iloc)), &
+                                                              &x_face_quad(:,ggi), normal(:,ggi))
+
+                           end if
+                           
+                           ! Form the face value = detwei * (relperm*absperm/visc)
+                           if (di%subcy_opt_sat%have .and. di%subcy_opt_sat%consistent) then
+
+                              face_value = detwei(ggi) * &
+                                           sum(di%cached_face_value%relperm(:,ggi,vele,p)) * &
+                                           absperm_ele(1) / &
+                                           visc_ele(1)
+                           
+                           else 
+                           
+                              face_value = detwei(ggi) * &
+                                           di%cached_face_value%relperm(1,ggi,vele,p) * &
+                                           absperm_ele(1) / &
+                                           visc_ele(1)
+                           
+                           end if
+                           
+                           ! Form the local matrix given by - n_i . sum_{phase} ( relperm*absperm/visc ) dP_1/dx_j
+                           do jloc = 1,di%pressure_mesh%shape%loc
+                              
+                              p_mat_local(iloc,jloc) = p_mat_local(iloc,jloc) - &
+                                                       sum(p_dshape(jloc, ggi, :)*normgi, 1) * &
+                                                       face_value
+
+                              p_mat_local(oloc,jloc) = p_mat_local(oloc,jloc) - &
+                                                       sum(p_dshape(jloc, ggi, :)*(-normgi), 1) * &
+                                                       face_value
+
+                           end do
+
+                           ! Add gravity term to rhs = - n_i . sum_{phase} ( relperm*absperm/visc ) * den*grav
+
+                           ! Find g dot n
+                           g_dot_n = dot_product(grav_ele(:,1), normgi)
+
+                           p_rhs_local(iloc)  = p_rhs_local(iloc) - &
+                                                face_value * &
+                                                di%cached_face_value%den(ggi,vele,p) * &
+                                                g_dot_n
+                           
+                           ! Add opposite node value with change of sign due to opposite normal direction
+                           p_rhs_local(oloc)  = p_rhs_local(oloc) + &
+                                                face_value * &
+                                                di%cached_face_value%den(ggi,vele,p) * &
+                                                g_dot_n
+
+                           ! Add capilliary pressure term to rhs = n_i . sum_{phase} ( relperm*absperm/visc ) dP_c/dx_j
+                           ! only for phase > 1
+                           if (p > 1) then
+
+                              ! Find grad_P_c dot n
+                              grad_cap_p_dot_n  = dot_product(grad_cap_pressure_face_quad(:,ggi), normgi)
+
+                              p_rhs_local(iloc) = p_rhs_local(iloc) + &
+                                                  face_value * &
+                                                  grad_cap_p_dot_n
+
+                              ! Add opposite node value with change of sign due to opposite normal direction
+                              p_rhs_local(oloc) = p_rhs_local(oloc) - &
+                                                  face_value * &
+                                                  grad_cap_p_dot_n
+
+                           end if
+
+                        end if check_visited
+
+                     end do quadrature_loop
+
+                  end if is_neigh
+
+               end do face_loop
+
+            end do nodal_loop_i
+            
+         end subroutine pressure_volume_element_local_assemble 
+
+      ! --------------------------------------------------------------------------------------
+      
+         subroutine pressure_surface_element_local_assemble(di) 
+
+            !!< Assemble the local surface element pressure matrix and rhs
+         
+            type(darcy_impes_type), intent(inout) :: di
+                        
+            p_matrix_local_bdy = 0.0
+            p_rhs_local_bdy    = 0.0
+
+            p_nodes_bdy = face_global_nodes(di%pressure_mesh, sele)
+            
+            ! A no_normal_flow BC adds nothing to the matrix and rhs so return zeros
+            if (di%v_bc_flag(sele) == V_BC_TYPE_NO_NORMAL_FLOW) return
+            
+            ! If phase 1 and strong pressure BC then no need to add integrals so return zeros
+            if (di%pressure_bc_flag(sele) == PRESSURE_BC_TYPE_DIRICHLET) return
+            
+            ! obtain the transformed determinant*weight and normals
+            if (di%cached_face_value%cached_detwei_normal) then
+
+               detwei_bdy => di%cached_face_value%detwei_bdy(:,sele)
+               
+               normal_bdy => di%cached_face_value%normal_bdy(:,:,sele)
+            
+            else
+
+               x_ele     = ele_val(di%positions, face_ele(di%positions, sele))
+               x_ele_bdy = face_val(di%positions, sele)
+            
+               call transform_cvsurf_facet_to_physical(x_ele, x_ele_bdy, di%x_cvbdyshape, normal_bdy, detwei_bdy)
+            
+            end if
+            
+            if (di%v_bc_flag(sele) == V_BC_TYPE_PRESCRIBED_NORMAL_FLOW) then
+                              
+               bc_sele_val = ele_val(di%v_bc_value, sele)
+                        
+            else
+            
+               ! get the viscosity value for this sele for this phase
+               visc_ele_bdy = face_val(di%viscosity(p)%ptr, sele)
+
+               ! get the absolute permeability value for this sele
+               absperm_ele_bdy = face_val(di%absolute_permeability, sele)         
+               
+               ! the inverse characteristic length used for pressure weak bc
+               inv_char_len_ele_bdy = ele_val(di%inverse_characteristic_length, sele)
+                              
+               ! Get the pressure BC value if required
+               if (di%pressure_bc_flag(sele) == PRESSURE_BC_TYPE_WEAKDIRICHLET) then
+               
+                  bc_sele_val = ele_val(di%pressure_bc_value, sele)
+               
+               end if
+
+               ! The gravity values for this element for each direction
+               grav_ele_bdy = face_val(di%gravity, sele) 
+
+!!!!! *** THIS IS NOT POSSIBLE YET ***
+!!!!!               ! get the gradient capilliary pressure at the cv surface quadrature points for each direction 
+!!!!!               if (p > 1) then
+!!!!!                  grad_cap_pressure_face_quad_bdy = &
+!!!!!                 &darcy_impes_face_grad_at_quad_scalar(di%capilliary_pressure(p)%ptr, sele, dn = p_dshape)
+!!!!!               end if
+
+!!!!!               ! obtain the derivative of the pressure mesh shape function at the CV face quadrature points
+!!!!!               if (di%cached_face_value%cached_p_dshape) then
+
+!!!!!                  p_dshape_bdy => di%cached_face_value%p_dshape_bdy(:,:,:,sele)
+
+!!!!!               else
+
+!!!!!                  call transform_to_physical(di%positions, sele, x_shape = di%x_cvbdyshape_full, &
+!!!!!                                             shape = di%p_cvbdyshape_full, dshape = p_dshape_bdy)
+
+!!!!!               end if
+               
+            end if 
+
+            bc_iloc_loop: do iloc = 1, di%pressure_mesh%faces%shape%loc
+
+               bc_face_loop: do face = 1, di%cvfaces%sfaces
+
+                  bc_neigh_if: if(di%cvfaces%sneiloc(iloc,face)/=0) then
+
+                     bc_quad_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                        ggi = (face-1)*di%cvfaces%shape%ngi + gi
+                                                
+                        prescribed_normal_flow: if (di%v_bc_flag(sele) == V_BC_TYPE_PRESCRIBED_NORMAL_FLOW) then
+                           
+                           ! If have prescribed_normal_flow BC for this phase then include CV integral of value
+                           
+                           if (di%subcy_opt_sat%have .and. di%subcy_opt_sat%consistent) then
+                                                      
+                              p_rhs_local_bdy(iloc) = p_rhs_local_bdy(iloc) - &
+                                                      di%subcy_opt_sat%number * &
+                                                      bc_sele_val(iloc) * &
+                                                      detwei_bdy(ggi)
+                           
+                           else 
+                           
+                              p_rhs_local_bdy(iloc) = p_rhs_local_bdy(iloc) - &
+                                                      bc_sele_val(iloc) * &
+                                                      detwei_bdy(ggi)
+                           
+                           end if
+                           
+                        else prescribed_normal_flow
+                                                                               
+                           ! Form the face value = detwei * (relperm*absperm/visc)
+                           if (di%subcy_opt_sat%have .and. di%subcy_opt_sat%consistent) then
+                           
+                              face_value = detwei_bdy(ggi) * &
+                                           sum(di%cached_face_value%relperm_bdy(:,ggi,sele,p)) * &
+                                           absperm_ele_bdy(1) / &
+                                           visc_ele_bdy(1) 
+                           
+                           else 
+                           
+                              face_value = detwei_bdy(ggi) * &
+                                           di%cached_face_value%relperm_bdy(1,ggi,sele,p) * &
+                                           absperm_ele_bdy(1) / &
+                                           visc_ele_bdy(1) 
+                           
+                           end if
+                           
+                           ! Add gravity term to rhs = - n_i . sum_{phase} ( relperm*absperm/visc ) * den*grav
+
+                           ! Find g dot n
+                           g_dot_n = dot_product(grav_ele_bdy(:,1), normal_bdy(:,ggi))
+                           
+                           p_rhs_local_bdy(iloc)  = p_rhs_local_bdy(iloc) - &
+                                                    face_value * &
+                                                    di%cached_face_value%den_bdy(ggi,sele,p) * &
+                                                    g_dot_n
+
+!!!!! *** THIS IS NOT POSSIBLE YET ***
+!!!!!                           ! Add capilliary pressure term to rhs = n_i . sum_{phase} ( relperm*absperm/visc ) dP_c/dx_j
+!!!!!                           ! only for phase > 1
+!!!!!                           if (p > 1) then
+
+!!!!!                              ! Find grad_P_c dot n
+!!!!!                              grad_cap_p_dot_n  = dot_product(grad_cap_pressure_face_quad_bdy(:,ggi), normal_bdy(:,ggi))
+
+!!!!!                              p_rhs_local(iloc) = p_rhs_local(iloc) + &
+!!!!!                                                  face_value * &
+!!!!!                                                  grad_cap_p_dot_n
+
+!!!!!                           end if
+                           
+!!!!!                           ! Add implicit gradient of pressure phase 1 to matrix with coeff 
+!!!!!                           ! - n_i . sum_{phase} ( relperm*absperm/visc ) 
+
+                           
+                           ! If have phase weak pressure BC then include integral of BC value in rhs
+                           ! and an implicit surface mass matrix term
+                           pressure_weak: if (di%pressure_bc_flag(sele) == PRESSURE_BC_TYPE_WEAKDIRICHLET) then
+                                                            
+                              do jloc = 1,di%pressure_mesh%faces%shape%loc 
+                              
+                                 p_rhs_local_bdy(iloc) = p_rhs_local_bdy(iloc) - &
+                                                         inv_char_len_ele_bdy(jloc) * &
+                                                         di%p_cvbdyshape%n(jloc,ggi) * &
+                                                         bc_sele_val(jloc) * &
+                                                         sum(normal_bdy(:,ggi)) * &
+                                                         face_value * &
+                                                         di%weak_pressure_bc_coeff
+                                                            
+                              end do
+                                                         
+                              do jloc = 1,di%pressure_mesh%faces%shape%loc 
+
+                                 p_matrix_local_bdy(iloc,jloc) = p_matrix_local_bdy(iloc,jloc) - &
+                                                                 inv_char_len_ele_bdy(jloc) * &
+                                                                 di%p_cvbdyshape%n(jloc,ggi) * &
+                                                                 sum(normal_bdy(:,ggi)) * &
+                                                                 face_value * &
+                                                                 di%weak_pressure_bc_coeff
+
+                              end do
+
+                           end if pressure_weak
+
+                        end if prescribed_normal_flow
+                        
+                     end do bc_quad_loop
+
+                  end if bc_neigh_if
+
+               end do bc_face_loop
+
+            end do bc_iloc_loop            
+            
+         end subroutine pressure_surface_element_local_assemble 
+
+      ! --------------------------------------------------------------------------------------
+      
+   end subroutine darcy_impes_assemble_and_solve_first_phase_pressure
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_get_v_boundary_condition(v, &
+                                                   types, &
+                                                   bc_surface_mesh, &
+                                                   v_bc_value, &
+                                                   v_bc_flag)
+
+      !!< Form the data associated with any v BC by returning 
+      !!< full surface mesh arrays of a flag indicating  either 
+      !!< no_normal_flow or prescribed_normal_flow and the value.
+
+      type(vector_field),               intent(in),   target :: v
+      character(len=*),   dimension(:), intent(in)           :: types
+      type(mesh_type),                  intent(in)           :: bc_surface_mesh
+      type(scalar_field),               intent(inout)        :: v_bc_value
+      integer,            dimension(:), intent(inout)        :: v_bc_flag
+
+      ! Local variables
+      character(len=FIELD_NAME_LEN)                       :: bctype
+      type(scalar_field)                                  :: scalar_surface_field
+      integer,                      dimension(:), pointer :: surface_element_list
+      integer                                             :: i, j, k, sele
+
+      ewrite(1,*) 'Get DarcyVelocity boundary condition data'
+
+      ! Zero the prescribed_normal_flow value surface field on whole boundary mesh  
+      call zero(v_bc_value)
+
+      ! Initialise flag for whether surface element has prescribed_normal_flow BC
+      v_bc_flag = 0
+
+      ! Loop each BC object instance for the v
+      ! (May have multiple prescribed_normal_flow BC's applied to different surface id's)
+      BC_loop: do i=1, get_boundary_condition_count(v)
+
+         ! Get this BC info
+         call get_boundary_condition(v, i, type = bctype, &
+                                     surface_element_list = surface_element_list)
+
+         ! check this is a prescribed_normal_flow BC  or no_normal_flow (nothing else is permitted)
+         if ((trim(bctype) /= 'prescribed_normal_flow') .and. (trim(bctype) /= 'no_normal_flow')) then
+            FLAbort('Have unknown BC type for a DarcyVelocity')
+         end if
+
+         ! Extract the surface_field for this BC for prescribed_normal_flow
+         if (trim(bctype) == 'prescribed_normal_flow') then
+            if (associated(v%bc%boundary_condition(i)%surface_fields)) then
+               scalar_surface_field = extract_scalar_field(&
+                    v%bc%boundary_condition(i)%surface_fields(1), 1)
+            else
+               FLAbort('Component surface_fields for DarcyVelocity BC type not associated')
+            end if
+         end if 
+         
+         ! Loop the surface elements associated with this BC instance
+         ! and place the required BC value in whole boundary field
+         ! for the prescribed_normal_flow BC type
+         BC_sele_loop: do k = 1, size(surface_element_list)
+
+            ! Find the whole domain surface element number
+            sele = surface_element_list(k)
+
+            ! Check that there is only 1 BC applied per surface element
+            if (v_bc_flag(sele) /= 0) then             
+               FLExit('Cannot apply more than 1 BC to a surface element for DarcyVelocity')
+            end if
+
+            ! Set the sele flag to indicate the BC type
+            do j = 1, size(types)
+               if (trim(types(j)) == trim(bctype)) exit
+            end do
+            if (j > size(types)) then
+               FLAbort('Cannot find no_normal_flow or prescribed_normal_flow bctype for DarcyVelocity')
+            end if
+            
+            v_bc_flag(sele) = j
+
+            ! Set the prescribed_normal_flow field values from this BC for its sele
+            if (trim(bctype) == 'prescribed_normal_flow') then
+               call set(v_bc_value, &
+                        ele_nodes(bc_surface_mesh, sele), &
+                        ele_val(scalar_surface_field, k))
+            end if
+            
+         end do BC_sele_loop
+
+      end do BC_loop
+
+      ewrite_minmax(v_bc_value)
+
+      ewrite(1,*) 'Finished get DarcyVelocity boundary condition data'
+
+   end subroutine darcy_impes_get_v_boundary_condition
+  
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_get_entire_sfield_boundary_condition(sfield, &
+                                                               types, &
+                                                               bc_surface_mesh, &
+                                                               sfield_bc_value, &
+                                                               sfield_bc_flag)
+      
+      !!< Get the entire BC data for either the saturation or pressure field
+      
+      type(scalar_field),               intent(in),   target :: sfield
+      character(len=*),   dimension(:), intent(in)           :: types
+      type(mesh_type),                  intent(in)           :: bc_surface_mesh
+      type(scalar_field),               intent(inout)        :: sfield_bc_value
+      integer,            dimension(:), intent(inout)        :: sfield_bc_flag
+      
+      ! local variables
+      character(len=FIELD_NAME_LEN)                       :: bctype
+      type(scalar_field),                         pointer :: surface_field
+      integer,                      dimension(:), pointer :: surface_element_list
+      integer                                             :: i, j, k, sele
+
+      ewrite(1,*) 'Get ',trim(sfield%name),' boundary condition data'
+      
+      ! zero the bc value field
+      call zero(sfield_bc_value)
+      
+      ! Initialise the bc type list for each sele
+      sfield_bc_flag = 0
+      
+      ! Loop each BC object for the sfield 
+      ! (May have multiple BC's applied to different suface id's)
+      BC_loop: do i=1, get_boundary_condition_count(sfield)
+         
+         ! Get this BC info
+         call get_boundary_condition(sfield, i, type = bctype, &
+                                     surface_element_list = surface_element_list)
+
+         ! see if we're interested in this one, if not skip it
+         do j = 1, size(types)
+            if (trim(types(j)) == trim(bctype)) exit
+         end do
+         if (j > size(types)) cycle
+         
+         ! Extract the surface field for this BC
+         if (associated(sfield%bc%boundary_condition(i)%surface_fields)) then
+            ! extract 1st surface field
+            surface_field => sfield%bc%boundary_condition(i)%surface_fields(1)
+         else
+            FLAbort('Component surface_fields for Saturation or Pressure BC type not associated')
+         end if
+
+         ! Loop the surface elements associated with this BC instance
+         ! and place the required BC value in whole boundary field
+         BC_sele_loop: do k = 1, size(surface_element_list)
+            
+            ! Find the whole domain surface element number
+            sele = surface_element_list(k)
+            
+            ! Check that there is only 1 BC applied per surface element
+            if (sfield_bc_flag(sele)/=0) then
+               ewrite(0,*) 'Requested types:', types
+               ewrite(0,*) 'Of these boundary condition types only one may be applied'
+               ewrite(0,*) 'to each surface element.'
+               ewrite(0,*) 'Surface element nr.:', sele
+               ewrite(0,*) 'has types', types(sfield_bc_flag(sele)), bctype
+               ewrite(0,*) 'on field: ', sfield%name
+               FLAbort('Issue with BC data for Saturation or Pressure')
+            end if
+            
+            ! Set the sele flag to indicate the BC type
+            sfield_bc_flag(sele) = j
+            
+            ! Set the BC sfield values from this BC for its sele
+            call set(sfield_bc_value, &
+                     ele_nodes(bc_surface_mesh, sele), &
+                     ele_val(surface_field, k))
+
+         end do BC_sele_loop
+         
+      end do BC_loop      
+
+      ewrite_minmax(sfield_bc_value)
+      
+      ewrite(1,*) 'Finished get ',trim(sfield%name),' boundary condition data'
+      
+   end subroutine darcy_impes_get_entire_sfield_boundary_condition
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_calculate_non_first_phase_pressures(di)
+      
+      !!< Calcuate the non first phase pressures which = first_phase_pressure + capilliary pressure
+      
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variables
+      integer :: p
+
+      ewrite(1,*) 'Calculate Non first phase Pressures'
+      
+      phase_loop: do p = 2, di%number_phase
+      
+         call set(di%pressure(p)%ptr, di%pressure(1)%ptr)
+         
+         call addto(di%pressure(p)%ptr, di%capilliary_pressure(p)%ptr)
+         
+         ewrite_minmax(di%pressure(p)%ptr) 
+         
+      end do phase_loop
+
+      ewrite(1,*) 'Finished Calculate Non first phase Pressures'
+
+   end subroutine darcy_impes_calculate_non_first_phase_pressures
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_calculate_average_pressure(di)
+      
+      !!< Calcuate the average phase pressure
+      
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variables
+      integer :: p
+
+      ewrite(1,*) 'Calculate AveragePressure'
+      
+      call zero(di%average_pressure)
+      
+      phase_loop: do p = 1, di%number_phase
+               
+         call addto(di%average_pressure, di%pressure(p)%ptr)
+                  
+      end do phase_loop
+      
+      call scale(di%average_pressure, 1.0/real(di%number_phase))
+      
+      ewrite_minmax(di%average_pressure) 
+
+      ewrite(1,*) 'Finished Calculate AveragePressure'
+
+   end subroutine darcy_impes_calculate_average_pressure
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_calculate_gradient_pressures(di)
+      
+      !!< Calculate the gradient pressures for each phase
+      
+      ! This routine requires optimisation and would become 
+      ! not required if ele_grad_at_quad is used for pressure 
+      ! and a face_grad_at_quad is made available
+
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variables
+      integer :: p
+      
+      ewrite(1,*) 'Calculate GradientPressure for each phase'
+      
+      phase_loop: do p = 1, di%number_phase
+      
+         call grad(di%pressure(p)%ptr, &
+                   di%positions, &
+                   di%gradient_pressure(p)%ptr)
+
+         ewrite_minmax(di%gradient_pressure(p)%ptr)
+      
+      end do phase_loop
+
+      ewrite(1,*) 'Finished Calculate GradientPressure for each phase'
+      
+   end subroutine darcy_impes_calculate_gradient_pressures
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_calculate_divergence_total_darcy_velocity(di)
+      
+      !!< Calculate the divergence of the total divergence velocity, 
+      !!< which may include terms summed over subcycles
+
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variables
+      integer :: vele, p, iloc, oloc, jloc, face, gi, ggi, sele, dim
+      real    :: darcy_vel_face_value_dot_n, v_over_relperm_dot_n, face_value
+      real,    dimension(1)            :: absperm_ele, visc_ele
+      real,    dimension(:,:), pointer :: grad_pressure_face_quad
+      real,    dimension(:,:), pointer :: grav_ele
+      real,    dimension(:,:), pointer :: x_ele
+      real,    dimension(:,:), pointer :: normal
+      real,    dimension(:),   pointer :: detwei
+      real,    dimension(:),   pointer :: normgi
+      logical, dimension(:),   pointer :: notvisited
+      real,    dimension(:),   pointer :: div_tvel_rhs_local
+      real,    dimension(:,:), pointer :: x_face_quad
+      integer, dimension(:),   pointer :: x_pmesh_nodes
+      integer, dimension(:),   pointer :: p_nodes      
+      real,    dimension(1)            :: absperm_ele_bdy, visc_ele_bdy
+      real,    dimension(:,:), pointer :: grav_ele_bdy
+      real,    dimension(:,:), pointer :: grad_pressure_face_quad_bdy
+      real,    dimension(:,:), pointer :: v_over_relperm_face_quad_bdy
+      real,    dimension(:),   pointer :: bc_sele_val
+      real,    dimension(:),   pointer :: pressure_ele_bdy
+      real,    dimension(:),   pointer :: inv_char_len_ele_bdy
+      real,    dimension(:,:), pointer :: normal_bdy
+      real,    dimension(:),   pointer :: detwei_bdy
+      real,    dimension(:,:), pointer :: x_ele_bdy
+      real,    dimension(:),   pointer :: div_tvel_rhs_local_bdy
+      integer, dimension(:),   pointer :: p_nodes_bdy
+      integer, parameter :: V_BC_TYPE_PRESCRIBED_NORMAL_FLOW = 1, V_BC_TYPE_NO_NORMAL_FLOW = 2
+      integer, parameter :: PRESSURE_BC_TYPE_WEAKDIRICHLET = 1
+                        
+      ewrite(1,*) 'Calculate the DivergenceTotalDarcyVelocity'
+      
+      ! allocate arrays used in assemble process - many assume
+      ! that all elements are the same type.
+      allocate(x_ele(di%ndim, ele_loc(di%positions,1)))      
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(normal(di%ndim,di%x_cvshape%ngi))
+         allocate(detwei(di%x_cvshape%ngi))      
+      end if
+      allocate(normgi(di%ndim))
+      allocate(notvisited(di%x_cvshape%ngi))
+      allocate(div_tvel_rhs_local(ele_loc(di%pressure_mesh,1)))
+      allocate(x_face_quad(di%ndim, di%x_cvshape%ngi))
+      allocate(grad_pressure_face_quad(di%ndim, di%p_cvshape%ngi))
+      allocate(grav_ele(di%ndim,1))
+      
+      allocate(grav_ele_bdy(di%ndim,1))
+      allocate(grad_pressure_face_quad_bdy(di%ndim, di%p_cvbdyshape%ngi))
+      allocate(v_over_relperm_face_quad_bdy(di%ndim, di%p_cvbdyshape%ngi))
+      allocate(bc_sele_val(face_loc(di%pressure_mesh,1)))
+      allocate(pressure_ele_bdy(face_loc(di%pressure_mesh,1)))
+      allocate(inv_char_len_ele_bdy(face_loc(di%pressure_mesh,1)))
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(detwei_bdy(di%x_cvbdyshape%ngi))
+         allocate(normal_bdy(di%ndim, di%x_cvbdyshape%ngi))
+      end if
+      allocate(div_tvel_rhs_local_bdy(face_loc(di%pressure_mesh,1)))      
+      allocate(x_ele_bdy(di%ndim, face_loc(di%positions,1)))      
+      allocate(p_nodes_bdy(face_loc(di%pressure_mesh,1)))
+
+      call zero(di%div_total_darcy_velocity)
+
+      ! Get the pressure BC - required if no v given and weak pressure dirichlet given for extra integrals
+      call darcy_impes_get_entire_sfield_boundary_condition(di%pressure(1)%ptr, &
+                                                            (/"weakdirichlet"/), &
+                                                            di%bc_surface_mesh, &
+                                                            di%pressure_bc_value, &
+                                                            di%pressure_bc_flag)
+
+      phase_loop: do p = 1, di%number_phase
+                       
+         vele_loop: do vele = 1, di%number_vele
+
+            ! The node indices of the pressure field
+            p_nodes => ele_nodes(di%pressure_mesh, vele)
+            
+            ! The gravity values for this element for each direction
+            grav_ele = ele_val(di%gravity, vele) 
+
+            ! get the viscosity value for this element for this phase
+            visc_ele = ele_val(di%viscosity(p)%ptr, vele)
+
+            ! get the absolute permeability value for this element
+            absperm_ele = ele_val(di%absolute_permeability, vele)         
+
+            ! get the latest gradient pressure at the cv surface quadrature points for each direction for this phase
+            grad_pressure_face_quad = ele_val_at_quad(di%gradient_pressure(p)%ptr, vele, di%gradp_cvshape)
+
+            ! obtain the transformed determinant*weight and normals
+            if (di%cached_face_value%cached_detwei_normal) then
+               
+               detwei => di%cached_face_value%detwei(:,vele)
+               
+               normal => di%cached_face_value%normal(:,:,vele)
+               
+            else
+
+               ! get the coordinate values for this element for each positions local node
+               x_ele = ele_val(di%positions, vele)         
+
+               ! get the coordinate values for this element for each quadrature point
+               x_face_quad = ele_val_at_quad(di%positions, vele, di%x_cvshape)
+
+               ! The node indices of the positions projected to the pressure mesh
+               x_pmesh_nodes => ele_nodes(di%positions_pressure_mesh, vele)
+            
+               call transform_cvsurf_to_physical(x_ele, di%x_cvshape, detwei, normal, di%cvfaces)
+            
+            end if
+
+            ! Initialise array for the quadrature points of this 
+            ! element for whether it has already been visited
+            notvisited = .true.
+
+            ! Initialise the local rhs to assemble for this element
+            div_tvel_rhs_local = 0.0
+
+            ! loop over local nodes within this element
+            nodal_loop_i: do iloc = 1, di%pressure_mesh%shape%loc
+
+               ! loop over CV faces internal to this element
+               face_loop: do face = 1, di%cvfaces%faces
+
+                  ! is this a face neighbouring iloc?
+                  is_neigh: if(di%cvfaces%neiloc(iloc, face) /= 0) then
+
+                     ! find the opposing local node across the CV face
+                     oloc = di%cvfaces%neiloc(iloc, face)
+
+                     ! loop over gauss points on face
+                     quadrature_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                        ! global gauss pt index for this element
+                        ggi = (face-1)*di%cvfaces%shape%ngi + gi
+
+                        ! check if this quadrature point has already been visited
+                        check_visited: if(notvisited(ggi)) then
+
+                           notvisited(ggi) = .false.
+
+                           if (di%cached_face_value%cached_detwei_normal) then
+
+                              ! cached normal has correct orientation already
+                              normgi = normal(:,ggi)
+
+                           else
+
+                              ! correct the orientation of the normal so it points away from iloc
+                              normgi = orientate_cvsurf_normgi(node_val(di%positions_pressure_mesh, x_pmesh_nodes(iloc)), &
+                                                              &x_face_quad(:,ggi), normal(:,ggi))
+
+                           end if
+
+                           ! Form the face value = detwei * (relperm*absperm/visc)
+                           if (di%subcy_opt_sat%have .and. di%subcy_opt_sat%consistent) then
+                           
+                              face_value = detwei(ggi) * &
+                                           sum(di%cached_face_value%relperm(:,ggi,vele,p)) * &
+                                           absperm_ele(1) / &
+                                           visc_ele(1) 
+                           
+                           else
+                           
+                              face_value = detwei(ggi) * &
+                                           di%cached_face_value%relperm(1,ggi,vele,p) * &
+                                           absperm_ele(1) / &
+                                           visc_ele(1) 
+                           
+                           end if
+                           
+                           darcy_vel_face_value_dot_n = - dot_product( normgi, &
+                          &face_value * (grad_pressure_face_quad(:,ggi) - di%cached_face_value%den(ggi,vele,p) * grav_ele(:,1)) )
+
+                           div_tvel_rhs_local(iloc) = div_tvel_rhs_local(iloc) + &
+                                                      darcy_vel_face_value_dot_n 
+
+                           div_tvel_rhs_local(oloc) = div_tvel_rhs_local(oloc) - &
+                                                      darcy_vel_face_value_dot_n 
+
+                        end if check_visited
+
+                     end do quadrature_loop
+
+                  end if is_neigh
+
+               end do face_loop
+
+            end do nodal_loop_i
+
+            call addto(di%div_total_darcy_velocity, p_nodes, div_tvel_rhs_local)
+
+         end do vele_loop
+
+         ! Get this phase v BC info - only for no_normal_flow and prescribed_normal_flow which is special as it is a scalar
+         call darcy_impes_get_v_boundary_condition(di%darcy_velocity(p)%ptr, &
+                                                   (/"prescribed_normal_flow", &
+                                                     "no_normal_flow        "/), &
+                                                   di%bc_surface_mesh, &
+                                                   di%v_bc_value, &
+                                                   di%v_bc_flag)
+                  
+         sele_loop: do sele = 1, di%number_sele
+            
+            if (di%v_bc_flag(sele) == V_BC_TYPE_NO_NORMAL_FLOW) cycle sele_loop
+
+            p_nodes_bdy = face_global_nodes(di%pressure_mesh, sele)
+
+            ! obtain the transformed determinant*weight and normals
+            if (di%cached_face_value%cached_detwei_normal) then
+
+               detwei_bdy => di%cached_face_value%detwei_bdy(:,sele)
+               
+               normal_bdy => di%cached_face_value%normal_bdy(:,:,sele)
+            
+            else
+
+               x_ele     = ele_val(di%positions, face_ele(di%positions, sele))
+               x_ele_bdy = face_val(di%positions, sele)
+            
+               call transform_cvsurf_facet_to_physical(x_ele, x_ele_bdy, di%x_cvbdyshape, normal_bdy, detwei_bdy)
+            
+            end if
+                     
+            if (di%v_bc_flag(sele) == V_BC_TYPE_PRESCRIBED_NORMAL_FLOW) then
+              
+               bc_sele_val = ele_val(di%v_bc_value, sele)
+            
+            else
+            
+               visc_ele_bdy         = face_val(di%viscosity(p)%ptr, sele)
+               absperm_ele_bdy      = face_val(di%absolute_permeability, sele)
+               inv_char_len_ele_bdy = ele_val(di%inverse_characteristic_length, sele)
+               if (di%pressure_bc_flag(sele) == PRESSURE_BC_TYPE_WEAKDIRICHLET) then
+                   bc_sele_val = ele_val(di%pressure_bc_value, sele)     
+               end if
+               grad_pressure_face_quad_bdy = face_val_at_quad(di%gradient_pressure(p)%ptr, sele, di%gradp_cvbdyshape)         
+               grav_ele_bdy                = face_val(di%gravity, sele) 
+               pressure_ele_bdy            = face_val(di%pressure(p)%ptr, sele) 
+
+               ! the latest DarcyVelocity/relperm at the quadrature points
+               do dim = 1,di%ndim
+
+                  v_over_relperm_face_quad_bdy(dim,:) =  &
+- (absperm_ele_bdy(1) / visc_ele_bdy(1)) * &
+(grad_pressure_face_quad_bdy(dim,:) - face_val_at_quad(di%density(p)%ptr, sele, di%p_cvbdyshape) * grav_ele_bdy(dim,1))   
+
+               end do
+
+            end if
+
+            ! Initialise the local rhs to assemble for this element
+            div_tvel_rhs_local_bdy = 0.0
+
+            bc_iloc_loop: do iloc = 1, di%pressure_mesh%faces%shape%loc
+
+               bc_face_loop: do face = 1, di%cvfaces%sfaces
+
+                  bc_neigh_if: if(di%cvfaces%sneiloc(iloc,face)/=0) then
+
+                     bc_quad_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                        ggi = (face-1)*di%cvfaces%shape%ngi + gi
+                                                
+                        prescribed_normal_flow: if (di%v_bc_flag(sele) == V_BC_TYPE_PRESCRIBED_NORMAL_FLOW) then
+
+                           if (di%subcy_opt_sat%have .and. di%subcy_opt_sat%consistent) then
+                           
+                              div_tvel_rhs_local_bdy(iloc) = div_tvel_rhs_local_bdy(iloc) + &
+                                                             di%subcy_opt_sat%number * &
+                                                             bc_sele_val(iloc) * &
+                                                             detwei_bdy(ggi)                   
+                           
+                           else
+                           
+                              div_tvel_rhs_local_bdy(iloc) = div_tvel_rhs_local_bdy(iloc) + &
+                                                             bc_sele_val(iloc) * &
+                                                             detwei_bdy(ggi)                   
+                           
+                           end if
+                           
+                        else prescribed_normal_flow
+                         
+                           ! NOTE this v_over_relperm_dot_n does not contain the relperm terms
+                           v_over_relperm_dot_n = dot_product(v_over_relperm_face_quad_bdy(:,ggi), normal_bdy(:,ggi))
+                           
+                           ! Adding this term includes the gradient of pressure, gradient of
+                           ! capilliary pressure and the gravity term. 
+                           if (di%subcy_opt_sat%have .and. di%subcy_opt_sat%consistent) then
+                           
+                              div_tvel_rhs_local_bdy(iloc) = div_tvel_rhs_local_bdy(iloc) + &
+                                                             sum(di%cached_face_value%relperm_bdy(:,ggi,sele,p)) * &
+                                                             v_over_relperm_dot_n * &
+                                                             detwei_bdy(ggi)
+                              
+                           else
+                              
+                              div_tvel_rhs_local_bdy(iloc) = div_tvel_rhs_local_bdy(iloc) + &
+                                                             di%cached_face_value%relperm_bdy(1,ggi,sele,p) * &
+                                                             v_over_relperm_dot_n * &
+                                                             detwei_bdy(ggi)
+                           
+                           end if
+                           
+                           ! Add two extra terms associated with weak pressure BC's
+                                                      
+                           ! Add two extra terms associated with weak pressure BC's
+                           pressure_weak: if (di%pressure_bc_flag(sele) == PRESSURE_BC_TYPE_WEAKDIRICHLET) then
+
+                              ! Form the face value = detwei * (relperm*absperm/visc)
+                              if (di%subcy_opt_sat%have .and. di%subcy_opt_sat%consistent) then
+                              
+                                 face_value = detwei_bdy(ggi) * &
+                                              sum(di%cached_face_value%relperm_bdy(:,ggi,sele,p)) * &
+                                              absperm_ele_bdy(1) / &
+                                              visc_ele_bdy(1) 
+                                 
+                              else
+                                 
+                                 face_value = detwei_bdy(ggi) * &
+                                              di%cached_face_value%relperm_bdy(1,ggi,sele,p) * &
+                                              absperm_ele_bdy(1) / &
+                                              visc_ele_bdy(1) 
+                              
+                              end if
+                              
+                              do jloc = 1, di%pressure_mesh%faces%shape%loc 
+
+                                 div_tvel_rhs_local_bdy(iloc) = div_tvel_rhs_local_bdy(iloc) - &
+                                                                inv_char_len_ele_bdy(jloc) * &
+                                                                di%p_cvbdyshape%n(jloc,ggi) * &
+                                                                bc_sele_val(jloc) * &
+                                                                sum(normal_bdy(:,ggi)) * &
+                                                                face_value * &
+                                                                di%weak_pressure_bc_coeff
+
+                              end do
+
+                              do jloc = 1, di%pressure_mesh%faces%shape%loc 
+
+                                 div_tvel_rhs_local_bdy(iloc) = div_tvel_rhs_local_bdy(iloc) + &
+                                                                inv_char_len_ele_bdy(jloc) * &
+                                                                di%p_cvbdyshape%n(jloc,ggi) * &
+                                                                pressure_ele_bdy(jloc) * &
+                                                                sum(normal_bdy(:,ggi)) * &
+                                                                face_value * &
+                                                                di%weak_pressure_bc_coeff
+
+                              end do
+
+                           end if pressure_weak
+                        
+                        end if prescribed_normal_flow
+                        
+                     end do bc_quad_loop
+
+                  end if bc_neigh_if
+
+               end do bc_face_loop
+
+            end do bc_iloc_loop
+
+            call addto(di%div_total_darcy_velocity, p_nodes_bdy, div_tvel_rhs_local_bdy)
+
+         end do sele_loop
+         
+      end do phase_loop
+      
+      call scale(di%div_total_darcy_velocity, di%inverse_cv_mass_pressure_mesh)
+      
+      ewrite_minmax(di%div_total_darcy_velocity)
+      
+      ! deallocate local variables as required
+      deallocate(x_ele)
+      if (di%cached_face_value%cached_detwei_normal) then
+         nullify(detwei)
+         nullify(normal)      
+      else
+         deallocate(detwei)
+         deallocate(normal)
+      end if
+      deallocate(normgi)
+      deallocate(notvisited)
+      deallocate(div_tvel_rhs_local)      
+      deallocate(x_face_quad)
+      deallocate(grad_pressure_face_quad)
+      deallocate(grav_ele)
+
+      deallocate(grav_ele_bdy)
+      deallocate(grad_pressure_face_quad_bdy)
+      deallocate(v_over_relperm_face_quad_bdy)
+      deallocate(bc_sele_val)
+      deallocate(pressure_ele_bdy)
+      deallocate(inv_char_len_ele_bdy)
+      if (di%cached_face_value%cached_detwei_normal) then
+         nullify(detwei_bdy)
+         nullify(normal_bdy)      
+      else
+         deallocate(detwei_bdy)
+         deallocate(normal_bdy)
+      end if
+      deallocate(div_tvel_rhs_local_bdy)
+      deallocate(x_ele_bdy)
+      deallocate(p_nodes_bdy)
+
+      ewrite(1,*) 'Finished Calculate the DivergenceTotalDarcyVelocity'
+            
+   end subroutine darcy_impes_calculate_divergence_total_darcy_velocity
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_assemble_and_solve_phase_saturations(di, &
+                                                               form_new_subcycle_relperm_face_values, &
+                                                               have_dual)
+      
+      !!< Assemble and solve the saturations for each subcycle.
+      !!< If form_new_subcycle_relperm_face_values is true
+      !!< and there are subcycles then new relperm face values are 
+      !!< formed. Phase 1 is special as it may be either solved 
+      !!< prognostically or calculated diagnostically from the 
+      !!< other phase saturation values that are solved first (for each subcycle).
+      
+      type(darcy_impes_type), intent(inout) :: di
+      logical,                intent(in)    :: form_new_subcycle_relperm_face_values
+      logical,                intent(in)    :: have_dual      
+
+      ! local variables
+      integer :: p, isub, isub_index
+      real    :: alpha_start, alpha_end
+
+      ewrite(1,*) 'Assemble and solve the saturations for each phase'
+                  
+      ! Solve the saturations for each subcycle via:            
+      !  - Form the relperm face values for each phase if required
+      !  - For each phase:
+      !     - Form the lhs = cv_mass_sfield_mesh_with_porosity / dt
+      !     - Form the rhs time = cv_mass_pressure_mesh_with_old_porosity * old_saturation_field / dt and add to rhs
+      !     - Add the s_source to rhs, which may include the dual perm transmissibility term
+      !     - Assemble and add the rhs adv to rhs 
+      !     - Apply any strong dirichlet BCs
+      !     - solve for latest saturations and copy to start subcycle saturations step
+      
+      ! Note: the porosity at the start and end of a subcycle time step
+      ! are linearly interpolated values from the main time step start and end
+
+      ! Deduce the subcycle time step size
+      if (di%subcy_opt_sat%have) then
+
+         di%dt_subcycle = di%dt/real(di%subcy_opt_sat%number)
+
+      else 
+
+         di%dt_subcycle = di%dt
+
+      end if 
+
+      ! Set the initial old saturation subcycle field
+      do p = 1, di%number_phase
+         call set(di%old_saturation_subcycle(p)%ptr, di%old_saturation(p)%ptr)
+      end do
+
+      ! Get the pressure BC - required if no v given and weak pressure dirichlet given for extra integrals
+      call darcy_impes_get_entire_sfield_boundary_condition(di%pressure(1)%ptr, &
+                                                            (/"weakdirichlet"/), &
+                                                            di%bc_surface_mesh, &
+                                                            di%pressure_bc_value, &
+                                                            di%pressure_bc_flag)
+
+      sub_loop: do isub = 1, di%subcy_opt_sat%number
+         
+         ewrite(1,*) 'Start subcycle ',isub
+         
+         ! form the start and end of subcycle dt
+         ! porosity linear interpolents
+         alpha_start = (isub - 1) / di%subcy_opt_sat%number
+         alpha_end   = isub / di%subcy_opt_sat%number
+         
+         ! Determine the index for cached face values depending on subcycling
+         if (di%subcy_opt_sat%have .and. di%subcy_opt_sat%consistent) then
+            isub_index = isub
+         else
+            isub_index = 1
+         end if            
+
+         ! If this is not the first subcycle then form the relperms 
+         ! for the start of this subcycle from the end of last subcycle sat
+         ! and face values if required
+         if (form_new_subcycle_relperm_face_values .and. (isub > 1)) then
+            call darcy_impes_calculate_relperm_fields(di)
+            
+            call darcy_impes_calculate_relperm_isub_face_values(di, isub_index)
+         end if
+       
+         phase_loop: do p = di%number_phase, 1, -1
+            
+            ewrite(1,*) 'Assemble and solve phase ',p
+            
+            ! If this is phase 1 and it is diagnostic then calculate it and exit loop
+            if ((p == 1) .and. di%phase_one_saturation_diagnostic) then
+               
+               call darcy_impes_calculate_phase_one_saturation_diagnostic(di)
+
+               ! Copy latest phase 1 solution to old subcycle field
+               call set(di%old_saturation_subcycle(1)%ptr, di%saturation(1)%ptr)
+               
+               exit phase_loop
+               
+            end if
+            
+            ! Allocate and get the BC data. If the BC is time dependent then
+            ! the end of overall time step values are used for all subcycles. 
+
+            ! Get this phase v BC info - only for no_normal_flow and prescribed_normal_flow which is special as it is a scalar
+            call darcy_impes_get_v_boundary_condition(di%darcy_velocity(p)%ptr, &
+                                                      (/"prescribed_normal_flow", &
+                                                        "no_normal_flow        "/), &
+                                                      di%bc_surface_mesh, &
+                                                      di%v_bc_value, &
+                                                      di%v_bc_flag)
+            
+            call zero(di%lhs)
+            call zero(di%rhs)
+            
+            call addto(di%lhs, di%cv_mass_pressure_mesh_with_porosity, scale = alpha_end / di%dt_subcycle)
+
+            call addto(di%lhs, di%cv_mass_pressure_mesh_with_old_porosity, scale = (1.0 - alpha_end) / di%dt_subcycle)
+
+            ! form the rhs time contribution and add
+            call addto(di%rhs, di%cv_mass_pressure_mesh_with_porosity, scale = alpha_start / di%dt_subcycle)
+
+            call addto(di%rhs, di%cv_mass_pressure_mesh_with_old_porosity, scale = (1.0 - alpha_start) / di%dt_subcycle)
+
+            call scale(di%rhs, di%old_saturation_subcycle(p)%ptr)
+
+            ! Add the saturation_source contribution
+            call compute_cv_mass(di%positions, di%cv_mass_pressure_mesh_with_source, di%saturation_source(p)%ptr)
+            
+            ! Should this include porosity ?!?!
+            call addto(di%rhs, di%cv_mass_pressure_mesh_with_source)
+            
+            ! Add the dual source term
+            if (have_dual) then               
+               call compute_cv_mass(di%positions, di%cv_mass_pressure_mesh_with_lambda_dual, di%transmissibility_lambda_dual(p)%ptr)
+
+               call set(di%rhs_dual, di%pressure_other_porous_media(p)%ptr)
+               call addto(di%rhs_dual, di%pressure(p)%ptr, scale = -1.0)            
+               call scale(di%rhs_dual, di%cv_mass_pressure_mesh_with_lambda_dual)
+               
+               call addto(di%rhs, di%rhs_dual)
+            end if
+            
+            ! assemble the rhs adv contribution and add
+            call darcy_impes_assemble_saturation_rhs_adv(di, &
+                                                         p, &
+                                                         isub_index)
+
+            ! apply strong dirichlet BC
+            call apply_dirichlet_conditions(di%lhs, di%rhs, di%saturation(p)%ptr)
+
+            ! Solve for the saturation
+            di%saturation(p)%ptr%val = di%rhs%val / di%lhs%val
+            
+            ! Update the halos
+            call halo_update(di%saturation(p)%ptr)
+            
+            ! Set the strong BC nodes to the values to be consistent
+            call set_dirichlet_consistent(di%saturation(p)%ptr) 
+
+            ! Copy latest solution to old subcycle
+            call set(di%old_saturation_subcycle(p)%ptr, di%saturation(p)%ptr)
+         
+         end do phase_loop
+         
+      end do sub_loop
+      
+      ewrite(1,*) 'Finished assemble and solve the saturations for each phase'
+            
+   end subroutine darcy_impes_assemble_and_solve_phase_saturations
+
+! ----------------------------------------------------------------------------
+                  
+   subroutine darcy_impes_assemble_saturation_rhs_adv(di, &
+                                                      p, &
+                                                      isub_index)
+
+      !!< Assemble the rhs advection contribtion for saturation phase p
+      
+      type(darcy_impes_type), intent(inout) :: di
+      integer,                intent(in)    :: p      
+      integer,                intent(in)    :: isub_index    
+
+      ! local variables
+      integer :: vele, iloc, oloc, jloc, face, gi, ggi, sele, dim
+      real    :: face_value, v_over_relperm_dot_n
+      real,    dimension(1)            :: absperm_ele, visc_ele
+      real,    dimension(:,:), pointer :: grad_pressure_face_quad
+      real,    dimension(:,:), pointer :: grav_ele
+      real,    dimension(:,:), pointer :: x_ele
+      real,    dimension(:,:), pointer :: normal
+      real,    dimension(:),   pointer :: detwei
+      real,    dimension(:),   pointer :: normgi
+      logical, dimension(:),   pointer :: notvisited
+      real,    dimension(:),   pointer :: s_rhs_local
+      real,    dimension(:,:), pointer :: x_face_quad
+      integer, dimension(:),   pointer :: x_pmesh_nodes
+      integer, dimension(:),   pointer :: p_nodes      
+      real,    dimension(1)            :: absperm_ele_bdy, visc_ele_bdy
+      real,    dimension(:,:), pointer :: grav_ele_bdy
+      real,    dimension(:,:), pointer :: grad_pressure_face_quad_bdy
+      real,    dimension(:,:), pointer :: v_over_relperm_face_quad_bdy
+      real,    dimension(:),   pointer :: bc_sele_val
+      real,    dimension(:),   pointer :: pressure_ele_bdy
+      real,    dimension(:),   pointer :: inv_char_len_ele_bdy
+      real,    dimension(:,:), pointer :: normal_bdy
+      real,    dimension(:),   pointer :: detwei_bdy
+      real,    dimension(:,:), pointer :: x_ele_bdy
+      real,    dimension(:),   pointer :: s_rhs_local_bdy
+      integer, dimension(:),   pointer :: p_nodes_bdy
+      integer, parameter :: V_BC_TYPE_PRESCRIBED_NORMAL_FLOW = 1, V_BC_TYPE_NO_NORMAL_FLOW = 2
+      integer, parameter :: PRESSURE_BC_TYPE_WEAKDIRICHLET = 1
+      
+      ! allocate arrays used in assemble process - many assume
+      ! that all elements are the same type.
+      allocate(x_ele(di%ndim, ele_loc(di%positions,1)))      
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(normal(di%ndim,di%x_cvshape%ngi))
+         allocate(detwei(di%x_cvshape%ngi))      
+      end if
+      allocate(normgi(di%ndim))
+      allocate(notvisited(di%x_cvshape%ngi))
+      allocate(s_rhs_local(ele_loc(di%pressure_mesh,1)))
+      allocate(x_face_quad(di%ndim, di%x_cvshape%ngi))
+      allocate(grad_pressure_face_quad(di%ndim, di%p_cvshape%ngi))
+      allocate(grav_ele(di%ndim,1))
+
+      allocate(grav_ele_bdy(di%ndim,1))
+      allocate(grad_pressure_face_quad_bdy(di%ndim, di%p_cvbdyshape%ngi))
+      allocate(v_over_relperm_face_quad_bdy(di%ndim, di%p_cvbdyshape%ngi))
+      allocate(bc_sele_val(face_loc(di%pressure_mesh,1)))
+      allocate(pressure_ele_bdy(face_loc(di%pressure_mesh,1)))
+      allocate(inv_char_len_ele_bdy(face_loc(di%pressure_mesh,1)))
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(detwei_bdy(di%x_cvbdyshape%ngi))
+         allocate(normal_bdy(di%ndim, di%x_cvbdyshape%ngi))
+      end if
+      allocate(s_rhs_local_bdy(face_loc(di%pressure_mesh,1)))
+      allocate(x_ele_bdy(di%ndim, face_loc(di%positions,1)))      
+      allocate(p_nodes_bdy(face_loc(di%pressure_mesh,1)))      
+      
+      ! Loop volume elements assembling local contributions    
+      vol_element_loop: do vele = 1, di%number_vele
+
+         ! The node indices of the pressure mesh (which saturation is associated with)
+         p_nodes => ele_nodes(di%pressure_mesh, vele)
+
+         ! The node indices of the positions projected to the pressure mesh (which saturation is associated with)
+         x_pmesh_nodes => ele_nodes(di%positions_pressure_mesh, vele)
+
+         ! The gravity values for this element for each direction
+         grav_ele = ele_val(di%gravity, vele) 
+
+         ! get the viscosity value for this element
+         visc_ele = ele_val(di%viscosity(p)%ptr, vele)
+
+         ! get the absolute permeability value for this element
+         absperm_ele = ele_val(di%absolute_permeability, vele)         
+
+         ! get the latest gradient pressure at the cv surface quadrature points for each direction 
+         grad_pressure_face_quad = ele_val_at_quad(di%gradient_pressure(p)%ptr, vele, di%gradp_cvshape)
+
+         ! obtain the transformed determinant*weight and normals
+         if (di%cached_face_value%cached_detwei_normal) then
+
+            detwei => di%cached_face_value%detwei(:,vele)
+
+            normal => di%cached_face_value%normal(:,:,vele)
+
+         else
+
+            ! get the coordinate values for this element for each positions local node
+            x_ele = ele_val(di%positions, vele)         
+
+            ! get the coordinate values for this element for each quadrature point
+            x_face_quad = ele_val_at_quad(di%positions, vele, di%x_cvshape)
+
+            call transform_cvsurf_to_physical(x_ele, di%x_cvshape, detwei, normal, di%cvfaces)
+
+         end if
+
+         ! Initialise array for the quadrature points of this 
+         ! element for whether it has already been visited
+         notvisited = .true.
+
+         ! Initialise the local rhs to assemble for this element
+         s_rhs_local = 0.0
+
+         ! loop over local nodes within this element
+         nodal_loop_i: do iloc = 1, di%pressure_mesh%shape%loc
+
+            ! loop over CV faces internal to this element
+            face_loop: do face = 1, di%cvfaces%faces
+
+               ! is this a face neighbouring iloc?
+               is_neigh: if(di%cvfaces%neiloc(iloc, face) /= 0) then
+
+                  ! find the opposing local node across the CV face
+                  oloc = di%cvfaces%neiloc(iloc, face)
+
+                  ! loop over gauss points on face
+                  quadrature_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                     ! global gauss pt index for this element
+                     ggi = (face-1)*di%cvfaces%shape%ngi + gi
+
+                     ! check if this quadrature point has already been visited
+                     check_visited: if(notvisited(ggi)) then
+
+                        notvisited(ggi) = .false.
+
+                        if (di%cached_face_value%cached_detwei_normal) then
+
+                           ! cached normal has correct orientation already
+                           normgi = normal(:,ggi)
+
+                        else
+
+                           ! correct the orientation of the normal so it points away from iloc
+                           normgi = orientate_cvsurf_normgi(node_val(di%positions_pressure_mesh, x_pmesh_nodes(iloc)), &
+                                                           &x_face_quad(:,ggi), normal(:,ggi))
+
+                        end if
+
+                        face_value = detwei(ggi) * di%cached_face_value%relperm(isub_index,ggi,vele,p) * absperm_ele(1) * &
+dot_product((grad_pressure_face_quad(:,ggi) - di%cached_face_value%den(ggi,vele,p) * grav_ele(:,1)), normgi)/ visc_ele(1)
+
+                        ! Form the local rhs for iloc and opposing oloc with normal vector sign change
+                        s_rhs_local(iloc) = s_rhs_local(iloc) + face_value
+
+                        s_rhs_local(oloc) = s_rhs_local(oloc) - face_value
+
+                     end if check_visited
+
+                  end do quadrature_loop
+
+               end if is_neigh
+
+            end do face_loop
+
+         end do nodal_loop_i
+
+         ! Add volume element contribution to global rhs field
+         call addto(di%rhs, p_nodes, s_rhs_local)
+
+      end do vol_element_loop
+
+      ! Add BC integrals
+
+      sele_loop: do sele = 1, di%number_sele
+
+         if (di%v_bc_flag(sele) == V_BC_TYPE_NO_NORMAL_FLOW) cycle sele_loop
+
+         p_nodes_bdy = face_global_nodes(di%pressure_mesh, sele)
+
+         ! obtain the transformed determinant*weight and normals
+         if (di%cached_face_value%cached_detwei_normal) then
+
+            detwei_bdy => di%cached_face_value%detwei_bdy(:,sele)
+
+            normal_bdy => di%cached_face_value%normal_bdy(:,:,sele)
+
+         else
+
+            x_ele     = ele_val(di%positions, face_ele(di%positions, sele))
+            x_ele_bdy = face_val(di%positions, sele)
+
+            call transform_cvsurf_facet_to_physical(x_ele, x_ele_bdy, di%x_cvbdyshape, normal_bdy, detwei_bdy)
+
+         end if
+
+         if (di%v_bc_flag(sele) == V_BC_TYPE_PRESCRIBED_NORMAL_FLOW) then
+
+            bc_sele_val = ele_val(di%v_bc_value, sele)
+
+         else
+
+            visc_ele_bdy         = face_val(di%viscosity(p)%ptr, sele)
+            absperm_ele_bdy      = face_val(di%absolute_permeability, sele)
+            inv_char_len_ele_bdy = ele_val(di%inverse_characteristic_length, sele)
+            if (di%pressure_bc_flag(sele) == PRESSURE_BC_TYPE_WEAKDIRICHLET) then
+                bc_sele_val = ele_val(di%pressure_bc_value, sele)     
+            end if
+            grad_pressure_face_quad_bdy = face_val_at_quad(di%gradient_pressure(p)%ptr, sele, di%gradp_cvbdyshape)         
+            grav_ele_bdy                = face_val(di%gravity, sele) 
+            pressure_ele_bdy            = face_val(di%pressure(p)%ptr, sele) 
+
+            ! the latest DarcyVelocity/relperm at the quadrature points for start of subcycle
+            ! determined from FE interpolation of each component, 
+            do dim = 1, di%ndim
+
+               v_over_relperm_face_quad_bdy(dim,:) =  &
+- (absperm_ele_bdy(1) / visc_ele_bdy(1)) * &
+(grad_pressure_face_quad_bdy(dim,:) - face_val_at_quad(di%density(p)%ptr, sele, di%p_cvbdyshape) * grav_ele_bdy(dim,1))   
+
+            end do
+
+         end if               
+
+         ! Initialise the local rhs to assemble for this element
+         s_rhs_local_bdy = 0.0
+
+         bc_iloc_loop: do iloc = 1, di%pressure_mesh%faces%shape%loc
+
+            bc_face_loop: do face = 1, di%cvfaces%sfaces
+
+               bc_neigh_if: if(di%cvfaces%sneiloc(iloc,face)/=0) then
+
+                  bc_quad_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                     ggi = (face-1)*di%cvfaces%shape%ngi + gi
+
+                     prescribed_normal_flow: if (di%v_bc_flag(sele) == V_BC_TYPE_PRESCRIBED_NORMAL_FLOW) then
+
+                        s_rhs_local_bdy(iloc) = s_rhs_local_bdy(iloc) - &
+                                                bc_sele_val(iloc) * &
+                                                detwei_bdy(ggi)
+
+                     else prescribed_normal_flow                           
+
+                        ! NOTE this v_over_relperm_dot_n does not contain the relperm terms
+                        v_over_relperm_dot_n = dot_product(v_over_relperm_face_quad_bdy(:,ggi), normal_bdy(:,ggi))
+
+                        ! This term includes the gradient pressure, gradient capilliary pressure 
+                        ! and the gravity term. 
+                        s_rhs_local_bdy(iloc) = s_rhs_local_bdy(iloc) - &
+                                                di%cached_face_value%relperm_bdy(isub_index,ggi,sele,p) * &
+                                                v_over_relperm_dot_n * &
+                                                detwei_bdy(ggi)
+
+                        ! Add two extra terms associated with weak pressure BC's                              
+                        pressure_weak: if (di%pressure_bc_flag(sele) == PRESSURE_BC_TYPE_WEAKDIRICHLET) then
+
+                           ! Form the face value = detwei * (relperm*absperm/visc)
+                           face_value = detwei_bdy(ggi) * &
+                                        di%cached_face_value%relperm_bdy(isub_index,ggi,sele,p) * &
+                                        absperm_ele_bdy(1) / &
+                                        visc_ele_bdy(1) 
+
+                           do jloc = 1, di%pressure_mesh%faces%shape%loc 
+
+                              s_rhs_local_bdy(iloc) = s_rhs_local_bdy(iloc) - &
+                                                      inv_char_len_ele_bdy(jloc) * &
+                                                      di%p_cvbdyshape%n(jloc,ggi) * &
+                                                      bc_sele_val(jloc) * &
+                                                      sum(normal_bdy(:,ggi)) * &
+                                                      face_value * &
+                                                      di%weak_pressure_bc_coeff
+
+                           end do
+
+                           do jloc = 1, di%pressure_mesh%faces%shape%loc 
+
+                              s_rhs_local_bdy(iloc) = s_rhs_local_bdy(iloc) + &
+                                                      inv_char_len_ele_bdy(jloc) * &
+                                                      di%p_cvbdyshape%n(jloc,ggi) * &
+                                                      pressure_ele_bdy(jloc) * &
+                                                      sum(normal_bdy(:,ggi)) * &
+                                                      face_value * &
+                                                      di%weak_pressure_bc_coeff
+
+                           end do
+
+                        end if pressure_weak
+
+                     end if prescribed_normal_flow
+
+                  end do bc_quad_loop
+
+               end if bc_neigh_if
+
+            end do bc_face_loop
+
+         end do bc_iloc_loop
+
+         call addto(di%rhs, p_nodes_bdy, s_rhs_local_bdy)
+
+      end do sele_loop
+
+      ! deallocate local variables as required
+      deallocate(x_ele)
+      if (di%cached_face_value%cached_detwei_normal) then
+         nullify(detwei)
+         nullify(normal)      
+      else
+         deallocate(detwei)
+         deallocate(normal)
+      end if
+      deallocate(normgi)
+      deallocate(notvisited)
+      deallocate(s_rhs_local)
+      deallocate(x_face_quad)
+      deallocate(grad_pressure_face_quad)
+      deallocate(grav_ele)
+
+      deallocate(grav_ele_bdy)
+      deallocate(grad_pressure_face_quad_bdy)
+      deallocate(v_over_relperm_face_quad_bdy)
+      deallocate(bc_sele_val)
+      deallocate(pressure_ele_bdy)
+      deallocate(inv_char_len_ele_bdy)
+      if (di%cached_face_value%cached_detwei_normal) then
+         nullify(detwei_bdy)
+         nullify(normal_bdy)      
+      else
+         deallocate(detwei_bdy)
+         deallocate(normal_bdy)
+      end if
+      deallocate(s_rhs_local_bdy)
+      deallocate(x_ele_bdy)
+      deallocate(p_nodes_bdy)
+
+   end subroutine darcy_impes_assemble_saturation_rhs_adv
+   
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_assemble_and_solve_generic_prog_sfields(di)
+      
+      !!< Assemble and solve generic prognostic scalar fields
+      
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variables
+      integer :: f, p
+      
+      ewrite(1,*) 'Assemble and solve generic prognostic scalar fields'
+      
+      sfield_loop: do f = 1, size(di%generic_prog_sfield)
+         
+         p = di%generic_prog_sfield(f)%phase
+         
+         call darcy_impes_assemble_and_solve_generic_prog_sfield(di, f, p)
+                  
+      end do sfield_loop
+
+      ewrite(1,*) 'Finished assemble and solve generic prognostic scalar fields'
+
+   end subroutine darcy_impes_assemble_and_solve_generic_prog_sfields
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_assemble_and_solve_generic_prog_sfield(di, f, p)
+      
+      !!< Assemble and a solve generic prognostic scalar field given by index f and p
+      
+      type(darcy_impes_type), intent(inout) :: di
+      integer,                intent(in)    :: f
+      integer,                intent(in)    :: p
+      
+      ! local variables
+      integer :: i
+            
+      ewrite(1,*) 'Assemble and solve sfield ',trim(di%generic_prog_sfield(f)%sfield%name),' of phase ',p
+
+      ! Get this phase v BC info - only for no_normal_flow and prescribed_normal_flow which is special as it is a scalar
+      call darcy_impes_get_v_boundary_condition(di%darcy_velocity(p)%ptr, &
+                                                (/"prescribed_normal_flow", &
+                                                  "no_normal_flow        "/), &
+                                                di%bc_surface_mesh, &
+                                                di%v_bc_value, &
+                                                di%v_bc_flag)
+
+      ! Get the pressure BC - required if no v given and weak pressure dirichlet given for extra integrals
+      call darcy_impes_get_entire_sfield_boundary_condition(di%pressure(1)%ptr, &
+                                                            (/"weakdirichlet"/), &
+                                                            di%bc_surface_mesh, &
+                                                            di%pressure_bc_value, &
+                                                            di%pressure_bc_flag)
+
+      ! Get the sfield BC
+      call darcy_impes_get_entire_sfield_boundary_condition(di%generic_prog_sfield(f)%sfield, &
+                                                            (/"weakdirichlet"/), &
+                                                            di%bc_surface_mesh, &
+                                                            di%sfield_bc_value, &
+                                                            di%sfield_bc_flag)
+            
+      call zero(di%lhs)         
+      call zero(di%matrix)
+      call zero(di%rhs)
+      call zero(di%rhs_full)
+      
+      ! Add porosity*saturation(absorption + 1/dt) to lhs 
+      if (di%generic_prog_sfield(f)%have_abs) then
+         
+         call addto(di%lhs, di%generic_prog_sfield(f)%sfield_abs)         
+         
+      end if
+      
+      call addto(di%lhs, 1.0/di%dt)            
+      call scale(di%lhs, di%cv_mass_pressure_mesh_with_porosity)
+      call scale(di%lhs, di%saturation(p)%ptr)
+            
+      ! Add old_porosity*old_saturation*old_sfield/dt to rhs      
+      call addto(di%rhs, 1.0/di%dt)
+      call scale(di%rhs, di%cv_mass_pressure_mesh_with_old_porosity)
+      call scale(di%rhs, di%old_saturation(p)%ptr)
+      call scale(di%rhs, di%generic_prog_sfield(f)%old_sfield)
+      
+      ! Add source to rhs   
+      if (di%generic_prog_sfield(f)%have_src) then
+         
+         call compute_cv_mass(di%positions, di%cv_mass_pressure_mesh_with_source, di%generic_prog_sfield(f)%sfield_src)
+         
+         call addto(di%rhs, di%cv_mass_pressure_mesh_with_source)
+         
+      end if
+            
+      ! Add diagonal lhs to matrix
+      call addto_diag(di%matrix, di%lhs)
+      
+      ! Add implicit low order advection and diffusion terms to matrix and rhs
+      if (di%generic_prog_sfield(f)%have_adv .or. di%generic_prog_sfield(f)%have_diff) then
+         
+         call darcy_impes_assemble_generic_prog_sfield_adv_diff(di, f, p)
+      
+      end if
+      
+      ! Solve for each face value iteration (default is 1)
+      do i = 1, di%generic_prog_sfield(f)%sfield_cv_options%number_face_value_iteration
+         
+         call set(di%rhs_full, di%rhs)
+         
+         ! Assmemble the high resolution rhs and to to rhs_full
+         if (di%generic_prog_sfield(f)%have_adv .and. &
+             di%generic_prog_sfield(f)%sfield_cv_options%facevalue == DARCY_IMPES_CV_FACEVALUE_FINITEELEMENT) then
+            
+            call zero(di%rhs_high_resolution)
+            
+            call darcy_impes_assemble_generic_prog_sfield_rhs_high_resolution(di, f, p)
+            
+            call addto(di%rhs_full, di%rhs_high_resolution)
+                        
+         end if
+         
+         ! Apply any strong dirichlet BC's
+         call apply_dirichlet_conditions(di%matrix, di%rhs_full, di%generic_prog_sfield(f)%sfield)
+
+         ! Solve the sfield
+         call petsc_solve(di%generic_prog_sfield(f)%sfield, di%matrix, di%rhs_full, di%state(p))
+
+         ! Set the strong BC nodes to the values to be consistent
+         call set_dirichlet_consistent(di%generic_prog_sfield(f)%sfield) 
+      
+      end do
+      
+      ewrite(1,*) 'Finished assemble and solve sfield ',trim(di%generic_prog_sfield(f)%sfield%name),' of phase ',p
+
+   end subroutine darcy_impes_assemble_and_solve_generic_prog_sfield
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_assemble_generic_prog_sfield_adv_diff(di, f, p)
+      
+      !!< Assemble the advection and diffusion terms for the 
+      !!< generic prognostic scalar field given by index f and p
+      
+      type(darcy_impes_type), intent(inout) :: di
+      integer,                intent(in)    :: f
+      integer,                intent(in)    :: p
+      
+      ! local variables
+      logical :: inflow
+      integer :: vele, iloc, oloc, jloc, face, gi, ggi, sele
+      real    :: face_value, v_dot_n, income, sat_face_value
+      real,    dimension(1)            :: absperm_ele, visc_ele, porosity_ele
+      real,    dimension(:,:),   pointer :: grad_pressure_face_quad
+      real,    dimension(:,:),   pointer :: grav_ele
+      real,    dimension(:,:,:), pointer :: diff_ele
+      real,    dimension(:),     pointer :: sat_ele            
+      real,    dimension(:,:),   pointer :: x_ele
+      real,    dimension(:,:,:), pointer :: p_dshape
+      real,    dimension(:,:),   pointer :: normal
+      real,    dimension(:),     pointer :: detwei
+      real,    dimension(:),     pointer :: normgi
+      logical, dimension(:),     pointer :: notvisited
+      real,    dimension(:,:),   pointer :: matrix_local
+      real,    dimension(:,:),   pointer :: x_face_quad
+      integer, dimension(:),     pointer :: x_pmesh_nodes
+      integer, dimension(:),     pointer :: p_nodes      
+      real,    dimension(1)              :: absperm_ele_bdy, visc_ele_bdy
+      real,    dimension(:,:),   pointer :: grav_ele_bdy
+      real,    dimension(:),     pointer :: ghost_sfield_ele_bdy      
+      real,    dimension(:,:),   pointer :: grad_pressure_face_quad_bdy
+      real,    dimension(:,:),   pointer :: v_over_relperm_face_quad_bdy
+      real,    dimension(:),     pointer :: bc_sele_val
+      real,    dimension(:),     pointer :: pressure_ele_bdy
+      real,    dimension(:),     pointer :: inv_char_len_ele_bdy
+      real,    dimension(:,:),   pointer :: normal_bdy
+      real,    dimension(:),     pointer :: detwei_bdy
+      real,    dimension(:,:),   pointer :: x_ele_bdy
+      real,    dimension(:),     pointer :: rhs_local_bdy
+      real,    dimension(:),     pointer :: matrix_local_bdy
+      integer, dimension(:),     pointer :: p_nodes_bdy
+      integer, parameter :: V_BC_TYPE_PRESCRIBED_NORMAL_FLOW = 1, V_BC_TYPE_NO_NORMAL_FLOW = 2
+      integer, parameter :: PRESSURE_BC_TYPE_WEAKDIRICHLET = 1
+      integer, parameter :: SFIELD_BC_TYPE_WEAKDIRICHLET   = 1
+      
+      ! allocate arrays used in assemble process - many assume
+      ! that all elements are the same type.
+      allocate(x_ele(di%ndim, ele_loc(di%positions,1)))      
+      if (.not. di%cached_face_value%cached_p_dshape) then
+         allocate(p_dshape(ele_loc(di%pressure_mesh,1), di%p_cvshape%ngi, mesh_dim(di%pressure_mesh)))
+      end if
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(normal(di%ndim,di%x_cvshape%ngi))
+         allocate(detwei(di%x_cvshape%ngi))      
+      end if
+      allocate(normgi(di%ndim))
+      allocate(notvisited(di%x_cvshape%ngi))
+      allocate(matrix_local(ele_loc(di%pressure_mesh,1),ele_loc(di%pressure_mesh,1)))
+      allocate(x_face_quad(di%ndim, di%x_cvshape%ngi))
+      allocate(grad_pressure_face_quad(di%ndim, di%p_cvshape%ngi))
+      allocate(grav_ele(di%ndim,1))
+      allocate(diff_ele(di%ndim,di%ndim,1))
+      allocate(sat_ele(ele_loc(di%positions,1)))      
+
+      allocate(grav_ele_bdy(di%ndim,1))
+      allocate(ghost_sfield_ele_bdy(face_loc(di%pressure_mesh,1)))
+      allocate(grad_pressure_face_quad_bdy(di%ndim, di%p_cvbdyshape%ngi))
+      allocate(v_over_relperm_face_quad_bdy(di%ndim, di%p_cvbdyshape%ngi))
+      allocate(bc_sele_val(face_loc(di%pressure_mesh,1)))
+      allocate(pressure_ele_bdy(face_loc(di%pressure_mesh,1)))
+      allocate(inv_char_len_ele_bdy(face_loc(di%pressure_mesh,1)))
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(detwei_bdy(di%x_cvbdyshape%ngi))
+         allocate(normal_bdy(di%ndim, di%x_cvbdyshape%ngi))
+      end if
+      allocate(rhs_local_bdy(face_loc(di%pressure_mesh,1)))
+      allocate(matrix_local_bdy(face_loc(di%pressure_mesh,1)))
+      allocate(x_ele_bdy(di%ndim, face_loc(di%positions,1)))      
+      allocate(p_nodes_bdy(face_loc(di%pressure_mesh,1)))      
+      
+      ! Loop volume elements assembling local contributions    
+      vol_element_loop: do vele = 1, di%number_vele
+
+         ! The node indices of the pressure mesh (which saturation is associated with)
+         p_nodes => ele_nodes(di%pressure_mesh, vele)
+
+         ! The node indices of the positions projected to the pressure mesh (which saturation is associated with)
+         x_pmesh_nodes => ele_nodes(di%positions_pressure_mesh, vele)
+         
+         ! The gravity values for this element for each direction
+         grav_ele = ele_val(di%gravity, vele) 
+
+         ! get the viscosity value for this element
+         visc_ele = ele_val(di%viscosity(p)%ptr, vele)
+
+         ! get the absolute permeability value for this element
+         absperm_ele = ele_val(di%absolute_permeability, vele)         
+
+         ! get the latest gradient pressure at the cv surface quadrature points for each direction 
+         grad_pressure_face_quad = ele_val_at_quad(di%gradient_pressure(p)%ptr, vele, di%gradp_cvshape)
+                          
+         ! obtain the transformed determinant*weight and normals
+         if (di%cached_face_value%cached_detwei_normal) then
+
+            detwei => di%cached_face_value%detwei(:,vele)
+
+            normal => di%cached_face_value%normal(:,:,vele)
+
+         else
+
+            ! get the coordinate values for this element for each positions local node
+            x_ele = ele_val(di%positions, vele)         
+
+            ! get the coordinate values for this element for each quadrature point
+            x_face_quad = ele_val_at_quad(di%positions, vele, di%x_cvshape)
+
+            call transform_cvsurf_to_physical(x_ele, di%x_cvshape, detwei, normal, di%cvfaces)
+
+         end if
+                  
+         if (di%generic_prog_sfield(f)%have_diff) then
+       
+            ! get the sfield diffusivity 
+            diff_ele = ele_val(di%generic_prog_sfield(f)%sfield_diff,vele)
+            
+            ! get the element values for porosity
+            porosity_ele = ele_val(di%porosity, vele)
+            
+            ! get the element values for saturation
+            sat_ele = ele_val(di%saturation(p)%ptr, vele)
+            
+            ! obtain the derivative of the pressure mesh shape function at the CV face quadrature points
+            if (di%cached_face_value%cached_p_dshape) then
+
+               p_dshape => di%cached_face_value%p_dshape(:,:,:,vele)
+
+            else
+
+               call transform_to_physical(di%positions, vele, x_shape = di%x_cvshape_full, &
+                                          shape = di%p_cvshape_full, dshape = p_dshape)
+
+            end if
+         
+         end if
+         
+         ! Initialise array for the quadrature points of this 
+         ! element for whether it has already been visited
+         notvisited = .true.
+
+         ! Initialise the local matrix to assemble for this element
+         matrix_local = 0.0
+
+         ! loop over local nodes within this element
+         nodal_loop_i: do iloc = 1, di%pressure_mesh%shape%loc
+
+            ! loop over CV faces internal to this element
+            face_loop: do face = 1, di%cvfaces%faces
+
+               ! is this a face neighbouring iloc?
+               is_neigh: if(di%cvfaces%neiloc(iloc, face) /= 0) then
+
+                  ! find the opposing local node across the CV face
+                  oloc = di%cvfaces%neiloc(iloc, face)
+
+                  ! loop over gauss points on face
+                  quadrature_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                     ! global gauss pt index for this element
+                     ggi = (face-1)*di%cvfaces%shape%ngi + gi
+
+                     ! check if this quadrature point has already been visited
+                     check_visited: if(notvisited(ggi)) then
+
+                        notvisited(ggi) = .false.
+
+                        if (di%cached_face_value%cached_detwei_normal) then
+
+                           ! cached normal has correct orientation already
+                           normgi = normal(:,ggi)
+
+                        else
+
+                           ! correct the orientation of the normal so it points away from iloc
+                           normgi = orientate_cvsurf_normgi(node_val(di%positions_pressure_mesh, x_pmesh_nodes(iloc)), &
+                                                           &x_face_quad(:,ggi), normal(:,ggi))
+
+                        end if
+
+                        v_dot_n = - di%cached_face_value%relperm(1,ggi,vele,p) * absperm_ele(1) * &
+dot_product((grad_pressure_face_quad(:,ggi) - di%cached_face_value%den(ggi,vele,p) * grav_ele(:,1)), normgi)/ visc_ele(1)
+
+                        inflow = (v_dot_n<=0.0)
+
+                        income = merge(1.0,0.0,inflow)                        
+                        
+                        if (di%generic_prog_sfield(f)%have_adv) then
+
+                           matrix_local(iloc, oloc) = matrix_local(iloc, oloc) + &
+                                                      detwei(ggi) * &
+                                                      v_dot_n * &
+                                                      income
+
+                           matrix_local(oloc, iloc) = matrix_local(oloc, iloc) + &
+                                                      detwei(ggi) * &
+                                                      (-v_dot_n) * &
+                                                      (1.0-income)
+
+                           matrix_local(iloc, iloc) = matrix_local(iloc, iloc) + &
+                                                      detwei(ggi) * &
+                                                      v_dot_n * &
+                                                      (1.0-income)
+
+                           matrix_local(oloc, oloc) = matrix_local(oloc, oloc) + &
+                                                      detwei(ggi) * &
+                                                      (-v_dot_n) * &
+                                                      income
+                        
+                        end if
+                        
+                        if (di%generic_prog_sfield(f)%have_diff) then
+                           
+                           ! Find the diffusion term saturation face value (taking upwind)
+                           sat_face_value = income*sat_ele(oloc) + (1.0-income)*sat_ele(iloc)
+                           
+                           do jloc=1, di%pressure_mesh%shape%loc
+                             
+                             matrix_local(iloc,jloc) = matrix_local(iloc,jloc) - &
+                                                       sum(matmul(diff_ele(:,:,1), p_dshape(jloc, ggi, :))*normgi, 1) * &
+                                                       porosity_ele(1) * &
+                                                       sat_face_value * &
+                                                       detwei(ggi)
+
+                             matrix_local(oloc, jloc) = matrix_local(oloc,jloc) - &
+                                                        sum(matmul(diff_ele(:,:,1), p_dshape(jloc, ggi, :))*(-normgi), 1) * &
+                                                        porosity_ele(1) * &
+                                                        sat_face_value * &
+                                                        detwei(ggi)
+                           
+                           end do                           
+                           
+                        end if
+                        
+                     end if check_visited
+
+                  end do quadrature_loop
+
+               end if is_neigh
+
+            end do face_loop
+
+         end do nodal_loop_i
+
+         call addto(di%matrix, p_nodes, p_nodes, matrix_local)
+
+      end do vol_element_loop
+
+      ! Add BC integrals associated with an advection term (diffusive flux is assumed zero on all boundary)
+      
+      have_adv_so_add_bc: if (di%generic_prog_sfield(f)%have_adv) then
+      
+         sele_loop: do sele = 1, di%number_sele
+
+            if (di%v_bc_flag(sele) == V_BC_TYPE_NO_NORMAL_FLOW) cycle sele_loop
+
+            p_nodes_bdy = face_global_nodes(di%pressure_mesh, sele)
+
+            ! obtain the transformed determinant*weight and normals
+            if (di%cached_face_value%cached_detwei_normal) then
+
+               detwei_bdy => di%cached_face_value%detwei_bdy(:,sele)
+
+               normal_bdy => di%cached_face_value%normal_bdy(:,:,sele)
+
+            else
+
+               x_ele     = ele_val(di%positions, face_ele(di%positions, sele))
+               x_ele_bdy = face_val(di%positions, sele)
+
+               call transform_cvsurf_facet_to_physical(x_ele, x_ele_bdy, di%x_cvbdyshape, normal_bdy, detwei_bdy)
+
+            end if
+
+            if (di%v_bc_flag(sele) == V_BC_TYPE_PRESCRIBED_NORMAL_FLOW) then
+
+               bc_sele_val = ele_val(di%v_bc_value, sele)
+
+            else
+
+               visc_ele_bdy         = face_val(di%viscosity(p)%ptr, sele)
+               absperm_ele_bdy      = face_val(di%absolute_permeability, sele)
+               inv_char_len_ele_bdy = ele_val(di%inverse_characteristic_length, sele)
+               if (di%pressure_bc_flag(sele) == PRESSURE_BC_TYPE_WEAKDIRICHLET) then
+                   bc_sele_val = ele_val(di%pressure_bc_value, sele)
+               end if
+               grad_pressure_face_quad_bdy = face_val_at_quad(di%gradient_pressure(p)%ptr, sele, di%gradp_cvbdyshape)         
+               grav_ele_bdy                = face_val(di%gravity, sele) 
+               pressure_ele_bdy            = face_val(di%pressure(p)%ptr, sele) 
+
+            end if
+
+            if (di%sfield_bc_flag(sele) == SFIELD_BC_TYPE_WEAKDIRICHLET) then
+               ghost_sfield_ele_bdy = ele_val(di%sfield_bc_value, sele)
+            else
+               ghost_sfield_ele_bdy = face_val(di%generic_prog_sfield(f)%sfield, sele)
+            end if
+
+            ! Initialise the local arrays to assemble for this element
+            rhs_local_bdy    = 0.0
+            matrix_local_bdy = 0.0
+
+            bc_iloc_loop: do iloc = 1, di%pressure_mesh%faces%shape%loc
+
+               bc_face_loop: do face = 1, di%cvfaces%sfaces
+
+                  bc_neigh_if: if(di%cvfaces%sneiloc(iloc,face)/=0) then
+
+                     bc_quad_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                        ggi = (face-1)*di%cvfaces%shape%ngi + gi
+
+                        prescribed_normal_flow: if (di%v_bc_flag(sele) == V_BC_TYPE_PRESCRIBED_NORMAL_FLOW) then
+
+                           if (bc_sele_val(iloc) <= 0) then
+
+                              income = 1.0
+
+                           else
+
+                              income = 0.0
+
+                           end if 
+
+                           rhs_local_bdy(iloc) = rhs_local_bdy(iloc) - &
+                                                 bc_sele_val(iloc) * &
+                                                 ghost_sfield_ele_bdy(iloc) * &
+                                                 detwei_bdy(ggi) * &
+                                                 income
+
+                           matrix_local_bdy(iloc) = matrix_local_bdy(iloc) + &
+                                                    bc_sele_val(iloc) * &
+                                                    detwei_bdy(ggi) * &
+                                                    (1.0 - income)
+
+                        else prescribed_normal_flow                          
+
+                           v_dot_n = - di%cached_face_value%relperm_bdy(1,ggi,sele,p) * absperm_ele_bdy(1) * &
+dot_product((grad_pressure_face_quad_bdy(:,ggi) - di%cached_face_value%den_bdy(ggi,sele,p) * grav_ele_bdy(:,1)), normal_bdy(:,ggi))/ visc_ele_bdy(1)
+
+                           if (v_dot_n <= 0) then
+
+                              income = 1.0
+
+                           else
+
+                              income = 0.0
+
+                           end if                         
+
+                           rhs_local_bdy(iloc) = rhs_local_bdy(iloc) - &
+                                                 v_dot_n * &
+                                                 ghost_sfield_ele_bdy(iloc) * &
+                                                 detwei_bdy(ggi) * &
+                                                 income
+
+                           matrix_local_bdy(iloc) = matrix_local_bdy(iloc) + &
+                                                    v_dot_n * &
+                                                    detwei_bdy(ggi) * &
+                                                    (1.0 - income)
+
+                           ! Add two extra terms associated with weak pressure BC's                              
+                           pressure_weak: if (di%pressure_bc_flag(sele) == PRESSURE_BC_TYPE_WEAKDIRICHLET) then
+
+                              ! Form the face value = detwei * (relperm*absperm/visc)
+                              face_value = detwei_bdy(ggi) * &
+                                           di%cached_face_value%relperm_bdy(1,ggi,sele,p) * &
+                                           absperm_ele_bdy(1) / &
+                                           visc_ele_bdy(1) 
+
+                              do jloc = 1, di%pressure_mesh%faces%shape%loc 
+
+                                 rhs_local_bdy(iloc) = rhs_local_bdy(iloc) - &
+                                                       inv_char_len_ele_bdy(jloc) * &
+                                                       di%p_cvbdyshape%n(jloc,ggi) * &
+                                                       bc_sele_val(jloc) * &
+                                                       sum(normal_bdy(:,ggi)) * &
+                                                       face_value * &
+                                                       ghost_sfield_ele_bdy(iloc) * &
+                                                       di%weak_pressure_bc_coeff * &
+                                                       income
+
+                                 matrix_local_bdy(iloc) = matrix_local_bdy(iloc) - &
+                                                          inv_char_len_ele_bdy(jloc) * &
+                                                          di%p_cvbdyshape%n(jloc,ggi) * &
+                                                          bc_sele_val(jloc) * &
+                                                          sum(normal_bdy(:,ggi)) * &
+                                                          face_value * &
+                                                          di%weak_pressure_bc_coeff * &
+                                                          (1.0 - income)
+
+                              end do
+
+                              do jloc = 1, di%pressure_mesh%faces%shape%loc 
+
+                                 rhs_local_bdy(iloc) = rhs_local_bdy(iloc) + &
+                                                       inv_char_len_ele_bdy(jloc) * &
+                                                       di%p_cvbdyshape%n(jloc,ggi) * &
+                                                       pressure_ele_bdy(jloc) * &
+                                                       sum(normal_bdy(:,ggi)) * &
+                                                       face_value * &
+                                                       ghost_sfield_ele_bdy(iloc) * &
+                                                       di%weak_pressure_bc_coeff * &
+                                                       income
+
+                                 matrix_local_bdy(iloc) = matrix_local_bdy(iloc) + &
+                                                          inv_char_len_ele_bdy(jloc) * &
+                                                          di%p_cvbdyshape%n(jloc,ggi) * &
+                                                          pressure_ele_bdy(jloc) * &
+                                                          sum(normal_bdy(:,ggi)) * &
+                                                          face_value * &
+                                                          di%weak_pressure_bc_coeff * &
+                                                          (1.0 - income)
+
+                              end do
+
+                           end if pressure_weak
+
+                        end if prescribed_normal_flow
+
+                     end do bc_quad_loop
+
+                  end if bc_neigh_if
+
+               end do bc_face_loop
+
+            end do bc_iloc_loop
+
+            call addto(di%rhs, p_nodes_bdy, rhs_local_bdy)
+
+            call addto_diag(di%matrix, p_nodes_bdy, matrix_local_bdy)
+
+         end do sele_loop
+      
+      end if have_adv_so_add_bc
+      
+      ! deallocate local variables as required
+      deallocate(x_ele)
+      if (di%cached_face_value%cached_p_dshape) then
+         nullify(p_dshape)      
+      else
+         deallocate(p_dshape)
+      end if
+      if (di%cached_face_value%cached_detwei_normal) then
+         nullify(detwei)
+         nullify(normal)      
+      else
+         deallocate(detwei)
+         deallocate(normal)
+      end if
+      deallocate(normgi)
+      deallocate(notvisited)
+      deallocate(matrix_local)
+      deallocate(x_face_quad)
+      deallocate(grad_pressure_face_quad)
+      deallocate(grav_ele)
+      deallocate(diff_ele)
+      deallocate(sat_ele)      
+
+      deallocate(grav_ele_bdy)
+      deallocate(ghost_sfield_ele_bdy)
+      deallocate(grad_pressure_face_quad_bdy)
+      deallocate(v_over_relperm_face_quad_bdy)
+      deallocate(bc_sele_val)
+      deallocate(pressure_ele_bdy)
+      deallocate(inv_char_len_ele_bdy)
+      if (di%cached_face_value%cached_detwei_normal) then
+         nullify(detwei_bdy)
+         nullify(normal_bdy)      
+      else
+         deallocate(detwei_bdy)
+         deallocate(normal_bdy)
+      end if
+      deallocate(rhs_local_bdy)
+      deallocate(matrix_local_bdy)
+      deallocate(x_ele_bdy)
+      deallocate(p_nodes_bdy)      
+      
+   end subroutine darcy_impes_assemble_generic_prog_sfield_adv_diff
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_assemble_generic_prog_sfield_rhs_high_resolution(di, f, p)
+      
+      !!< Assemble the advection high resolution terms for the 
+      !!< generic prognostic scalar field given by index f and p
+      
+      type(darcy_impes_type), intent(inout) :: di
+      integer,                intent(in)    :: f
+      integer,                intent(in)    :: p
+      
+      ! local variables
+      logical :: inflow
+      integer :: vele, iloc, oloc, face, gi, ggi, s_upwind_pos
+      real    :: v_dot_n, income, sfield_face_value, sfield_low_face_value, sfield_high_face_value
+      real,    dimension(1)            :: absperm_ele, visc_ele
+      real,    dimension(:,:),   pointer :: grad_pressure_face_quad
+      real,    dimension(:,:),   pointer :: grav_ele
+      real,    dimension(:),     pointer :: sfield_ele            
+      real,    dimension(:,:),   pointer :: x_ele
+      real,    dimension(:,:,:), pointer :: p_dshape
+      real,    dimension(:,:),   pointer :: normal
+      real,    dimension(:),     pointer :: detwei
+      real,    dimension(:),     pointer :: normgi
+      logical, dimension(:),     pointer :: notvisited
+      real,    dimension(:),     pointer :: rhs_local
+      real,    dimension(:,:),   pointer :: x_face_quad
+      integer, dimension(:),     pointer :: x_pmesh_nodes
+      integer, dimension(:),     pointer :: p_nodes      
+      integer, dimension(:),     pointer :: s_upwind_nodes
+      
+      ! allocate arrays used in assemble process - many assume
+      ! that all elements are the same type.
+      allocate(x_ele(di%ndim, ele_loc(di%positions,1)))      
+      if (.not. di%cached_face_value%cached_p_dshape) then
+         allocate(p_dshape(ele_loc(di%pressure_mesh,1), di%p_cvshape%ngi, mesh_dim(di%pressure_mesh)))
+      end if
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(normal(di%ndim,di%x_cvshape%ngi))
+         allocate(detwei(di%x_cvshape%ngi))      
+      end if
+      allocate(normgi(di%ndim))
+      allocate(notvisited(di%x_cvshape%ngi))
+      allocate(rhs_local(ele_loc(di%pressure_mesh,1)))
+      allocate(x_face_quad(di%ndim, di%x_cvshape%ngi))
+      allocate(grad_pressure_face_quad(di%ndim, di%p_cvshape%ngi))
+      allocate(grav_ele(di%ndim,1))
+      allocate(sfield_ele(ele_loc(di%positions,1)))      
+      
+      if(di%generic_prog_sfield(f)%sfield_cv_options%limit_facevalue) then
+
+        ! NOTE only new upwind values are required but this procedure 
+        !      expects latest and old. So pass in new twice - not optimal
+        call find_upwind_values(di%state, &
+                                di%positions_pressure_mesh, &
+                                di%generic_prog_sfield(f)%sfield, &
+                                di%sfield_upwind, &
+                                di%generic_prog_sfield(f)%sfield, &
+                                di%sfield_upwind, &
+                                option_path = trim(di%generic_prog_sfield(f)%sfield%option_path))
+
+      end if
+      
+      ! Initialise optimisation flag used in finding upwind value in high resolution schemes
+      s_upwind_pos = 0
+      
+      ! Loop volume elements assembling local contributions    
+      vol_element_loop: do vele = 1, di%number_vele
+
+         ! The node indices of the pressure mesh (which saturation is associated with)
+         p_nodes => ele_nodes(di%pressure_mesh, vele)
+
+         ! The node indices of the positions projected to the pressure mesh (which saturation is associated with)
+         x_pmesh_nodes => ele_nodes(di%positions_pressure_mesh, vele)
+         
+         ! The gravity values for this element for each direction
+         grav_ele = ele_val(di%gravity, vele) 
+
+         ! get the viscosity value for this element
+         visc_ele = ele_val(di%viscosity(p)%ptr, vele)
+
+         ! get the absolute permeability value for this element
+         absperm_ele = ele_val(di%absolute_permeability, vele)         
+
+         ! get the latest gradient pressure at the cv surface quadrature points for each direction 
+         grad_pressure_face_quad = ele_val_at_quad(di%gradient_pressure(p)%ptr, vele, di%gradp_cvshape)
+
+         ! get the sfield value for this element
+         sfield_ele = ele_val(di%generic_prog_sfield(f)%sfield, vele)
+
+         ! Determine the node numbers to use to determine the sfield upwind values
+         if((di%generic_prog_sfield(f)%sfield_cv_options%upwind_scheme == CV_UPWINDVALUE_PROJECT_POINT).or.&
+            (di%generic_prog_sfield(f)%sfield_cv_options%upwind_scheme == CV_UPWINDVALUE_PROJECT_GRAD)) then
+
+            s_upwind_nodes => x_pmesh_nodes
+
+         else
+
+            s_upwind_nodes => p_nodes
+
+         end if
+                          
+         ! obtain the transformed determinant*weight and normals
+         if (di%cached_face_value%cached_detwei_normal) then
+
+            detwei => di%cached_face_value%detwei(:,vele)
+
+            normal => di%cached_face_value%normal(:,:,vele)
+
+         else
+
+            ! get the coordinate values for this element for each positions local node
+            x_ele = ele_val(di%positions, vele)         
+
+            ! get the coordinate values for this element for each quadrature point
+            x_face_quad = ele_val_at_quad(di%positions, vele, di%x_cvshape)
+
+            call transform_cvsurf_to_physical(x_ele, di%x_cvshape, detwei, normal, di%cvfaces)
+
+         end if
+                  
+         ! Initialise array for the quadrature points of this 
+         ! element for whether it has already been visited
+         notvisited = .true.
+
+         ! Initialise the local rhs to assemble for this element
+         rhs_local = 0.0
+
+         ! loop over local nodes within this element
+         nodal_loop_i: do iloc = 1, di%pressure_mesh%shape%loc
+
+            ! loop over CV faces internal to this element
+            face_loop: do face = 1, di%cvfaces%faces
+
+               ! is this a face neighbouring iloc?
+               is_neigh: if(di%cvfaces%neiloc(iloc, face) /= 0) then
+
+                  ! find the opposing local node across the CV face
+                  oloc = di%cvfaces%neiloc(iloc, face)
+
+                  ! loop over gauss points on face
+                  quadrature_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                     ! global gauss pt index for this element
+                     ggi = (face-1)*di%cvfaces%shape%ngi + gi
+
+                     ! check if this quadrature point has already been visited
+                     check_visited: if(notvisited(ggi)) then
+
+                        notvisited(ggi) = .false.
+
+                        if (di%cached_face_value%cached_detwei_normal) then
+
+                           ! cached normal has correct orientation already
+                           normgi = normal(:,ggi)
+
+                        else
+
+                           ! correct the orientation of the normal so it points away from iloc
+                           normgi = orientate_cvsurf_normgi(node_val(di%positions_pressure_mesh, x_pmesh_nodes(iloc)), &
+                                                           &x_face_quad(:,ggi), normal(:,ggi))
+
+                        end if
+
+                        v_dot_n = - di%cached_face_value%relperm(1,ggi,vele,p) * absperm_ele(1) * &
+dot_product((grad_pressure_face_quad(:,ggi) - di%cached_face_value%den(ggi,vele,p) * grav_ele(:,1)), normgi)/ visc_ele(1)
+
+                        inflow = (v_dot_n<=0.0)
+
+                        income = merge(1.0,0.0,inflow)                        
+
+                        ! form the high resolution face value
+                        call darcy_impes_evaluate_sfield_face_val(sfield_high_face_value, &
+                                                                  iloc, &
+                                                                  oloc, &
+                                                                  ggi, &
+                                                                  s_upwind_nodes, &
+                                                                  di%p_cvshape, &
+                                                                  sfield_ele, &
+                                                                  di%sfield_upwind, &
+                                                                  inflow, &
+                                                                  income, &
+                                                                  di%generic_prog_sfield(f)%sfield_cv_options, &
+                                                                  save_pos = s_upwind_pos)                        
+                                                
+                        ! form and add the lower order iterated term to rhs
+                        sfield_low_face_value = income*sfield_ele(oloc) + (1.0-income)*sfield_ele(iloc)
+                        
+                        ! form the rhs face value and add to the rhs field
+                        sfield_face_value = sfield_low_face_value - sfield_high_face_value
+                        
+                        rhs_local(iloc) = rhs_local(iloc) + &
+                                          sfield_face_value * &
+                                          v_dot_n * &
+                                          detwei(ggi)
+
+                        rhs_local(oloc) = rhs_local(oloc)  - &
+                                          sfield_face_value * &
+                                          v_dot_n * &
+                                          detwei(ggi)
+                                                                       
+                     end if check_visited
+
+                  end do quadrature_loop
+
+               end if is_neigh
+
+            end do face_loop
+
+         end do nodal_loop_i
+
+         call addto(di%rhs_high_resolution, p_nodes, rhs_local)
+
+      end do vol_element_loop
+      
+      ! deallocate local variables as required
+      deallocate(x_ele)
+      if (di%cached_face_value%cached_p_dshape) then
+         nullify(p_dshape)      
+      else
+         deallocate(p_dshape)
+      end if
+      if (di%cached_face_value%cached_detwei_normal) then
+         nullify(detwei)
+         nullify(normal)      
+      else
+         deallocate(detwei)
+         deallocate(normal)
+      end if
+      deallocate(normgi)
+      deallocate(notvisited)
+      deallocate(rhs_local)
+      deallocate(x_face_quad)
+      deallocate(grad_pressure_face_quad)
+      deallocate(grav_ele)
+      deallocate(sfield_ele)      
+      
+   end subroutine darcy_impes_assemble_generic_prog_sfield_rhs_high_resolution
+
+! ----------------------------------------------------------------------------
+   
+   subroutine darcy_impes_calculate_phase_one_saturation_diagnostic(di)
+      
+      !!< Calculate the first phase diagnostic saturation as 1.0 - the rest
+      
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variables
+      integer :: p
+      
+      ewrite(1,*) 'Calculate phase 1 diagnostic Saturation'
+      
+      call set(di%saturation(1)%ptr, 1.0)
+      
+      do p = 2, di%number_phase
+         
+         call addto(di%saturation(1)%ptr, di%saturation(p)%ptr, scale = -1.0)
+         
+      end do
+      
+      ewrite_minmax(di%saturation(1)%ptr)
+       
+   end subroutine darcy_impes_calculate_phase_one_saturation_diagnostic
+
+! ----------------------------------------------------------------------------
+   
+   subroutine darcy_impes_calculate_vel_mob_ff_and_cfl_fields(di)
+      
+      !!< Calculate the various velocities, mobilities, fractional flows, and CFL fields
+      
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variables
+      type(vector_field) :: aux_vfield
+      logical :: inflow
+      integer :: vele, p, dim, iloc, oloc, face, gi, ggi, sele
+      real    :: income, v_dot_n
+      real,    dimension(1)              :: visc_ele, absperm_ele
+      real,    dimension(:,:),   pointer :: grad_pressure_face_quad
+      real,    dimension(:,:),   pointer :: v_local
+      real,    dimension(:),     pointer :: m_local
+      real,    dimension(:,:),   pointer :: grav_ele
+      real,    dimension(:,:),   pointer :: x_ele
+      real,    dimension(:,:),   pointer :: normal
+      real,    dimension(:),     pointer :: detwei
+      real,    dimension(:),     pointer :: normgi
+      logical, dimension(:),     pointer :: notvisited
+      real,    dimension(:),     pointer :: cfl_rhs_local
+      real,    dimension(:,:),   pointer :: x_face_quad
+      integer, dimension(:),     pointer :: x_pmesh_nodes
+      integer, dimension(:),     pointer :: p_nodes      
+      real,    dimension(1)              :: visc_ele_bdy, absperm_ele_bdy
+      real,    dimension(:,:),   pointer :: grav_ele_bdy
+      real,    dimension(:,:),   pointer :: grad_pressure_face_quad_bdy
+      real,    dimension(:,:),   pointer :: normal_bdy
+      real,    dimension(:),     pointer :: detwei_bdy
+      real,    dimension(:,:),   pointer :: x_ele_bdy
+      real,    dimension(:),     pointer :: cfl_rhs_local_bdy
+      real,    dimension(:,:),   pointer :: v_local_bdy
+      real,    dimension(:),     pointer :: m_local_bdy
+      integer, dimension(:),     pointer :: p_nodes_bdy
+      real,    dimension(:),     pointer :: bc_sele_val
+      integer, parameter :: V_BC_TYPE_PRESCRIBED_NORMAL_FLOW = 1, V_BC_TYPE_NO_NORMAL_FLOW = 2
+            
+      ewrite(1,*) 'Calculate CFL, Velocity, Mobilities and FractionalFlow fields'
+
+      ! allocate arrays used in CFL assemble process - many assume
+      ! that all elements are the same type.
+      allocate(x_ele(di%ndim, ele_loc(di%positions,1)))      
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(normal(di%ndim,di%x_cvshape%ngi))
+         allocate(detwei(di%x_cvshape%ngi))      
+      end if
+      allocate(normgi(di%ndim))
+      allocate(notvisited(di%x_cvshape%ngi))
+      allocate(cfl_rhs_local(ele_loc(di%pressure_mesh,1)))
+      allocate(x_face_quad(di%ndim, di%x_cvshape%ngi))
+      allocate(grad_pressure_face_quad(di%ndim, di%p_cvshape%ngi))
+      allocate(grav_ele(di%ndim,1))
+      
+      allocate(grav_ele_bdy(di%ndim,1))
+      allocate(grad_pressure_face_quad_bdy(di%ndim, di%p_cvbdyshape%ngi))
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(detwei_bdy(di%x_cvbdyshape%ngi))
+         allocate(normal_bdy(di%ndim, di%x_cvbdyshape%ngi))
+      end if
+      allocate(cfl_rhs_local_bdy(face_loc(di%pressure_mesh,1)))
+      allocate(x_ele_bdy(di%ndim, face_loc(di%positions,1)))      
+      allocate(p_nodes_bdy(face_loc(di%pressure_mesh,1)))
+      allocate(bc_sele_val(face_loc(di%pressure_mesh,1)))
+      
+      allocate(v_local(di%ndim,ele_loc(di%pressure_mesh,1)))
+      allocate(m_local(ele_loc(di%pressure_mesh,1)))
+      allocate(v_local_bdy(di%ndim,face_loc(di%pressure_mesh,1)))
+      allocate(m_local_bdy(face_loc(di%pressure_mesh,1)))
+
+      phase_loop: do p = 1, di%number_phase
+         
+         call zero(di%cfl(p)%ptr)
+
+         call zero(di%darcy_velocity(p)%ptr)
+         
+         call zero(di%mobility(p)%ptr)
+                               
+         vele_loop: do vele = 1, di%number_vele
+
+            ! The node indices of the pressure field
+            p_nodes => ele_nodes(di%pressure_mesh, vele)
+           
+            ! The gravity values for this element for each direction
+            grav_ele = ele_val(di%gravity, vele) 
+
+            ! get the viscosity value for this element for this phase
+            visc_ele = ele_val(di%viscosity(p)%ptr, vele)
+
+            ! get the absolute permeability value for this element
+            absperm_ele = ele_val(di%absolute_permeability, vele)         
+
+            ! get the latest gradient pressure at the cv surface quadrature points for each direction for this phase
+            grad_pressure_face_quad = ele_val_at_quad(di%gradient_pressure(p)%ptr, vele, di%gradp_cvshape)
+
+            ! obtain the transformed determinant*weight and normals
+            if (di%cached_face_value%cached_detwei_normal) then
+               
+               detwei => di%cached_face_value%detwei(:,vele)
+               
+               normal => di%cached_face_value%normal(:,:,vele)
+               
+            else
+
+               ! get the coordinate values for this element for each positions local node
+               x_ele = ele_val(di%positions, vele)         
+
+               ! get the coordinate values for this element for each quadrature point
+               x_face_quad = ele_val_at_quad(di%positions, vele, di%x_cvshape)
+
+               ! The node indices of the positions projected to the pressure mesh
+               x_pmesh_nodes => ele_nodes(di%positions_pressure_mesh, vele)
+            
+               call transform_cvsurf_to_physical(x_ele, di%x_cvshape, detwei, normal, di%cvfaces)
+            
+            end if
+
+            ! Initialise array for the quadrature points of this 
+            ! element for whether it has already been visited
+            notvisited = .true.
+
+            ! Initialise the local assemble arrays for this element
+            cfl_rhs_local = 0.0
+            v_local       = 0.0
+            m_local       = 0.0
+                        
+            ! loop over local nodes within this element
+            nodal_loop_i: do iloc = 1, di%pressure_mesh%shape%loc
+
+               ! loop over CV faces internal to this element
+               face_loop: do face = 1, di%cvfaces%faces
+
+                  ! is this a face neighbouring iloc?
+                  is_neigh: if(di%cvfaces%neiloc(iloc, face) /= 0) then
+
+                     ! find the opposing local node across the CV face
+                     oloc = di%cvfaces%neiloc(iloc, face)
+
+                     ! loop over gauss points on face
+                     quadrature_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                        ! global gauss pt index for this element
+                        ggi = (face-1)*di%cvfaces%shape%ngi + gi
+
+                        ! check if this quadrature point has already been visited
+                        check_visited: if(notvisited(ggi)) then
+
+                           notvisited(ggi) = .false.
+
+                           if (di%cached_face_value%cached_detwei_normal) then
+
+                              ! cached normal has correct orientation already
+                              normgi = normal(:,ggi)
+
+                           else
+
+                              ! correct the orientation of the normal so it points away from iloc
+                              normgi = orientate_cvsurf_normgi(node_val(di%positions_pressure_mesh, x_pmesh_nodes(iloc)), &
+                                                              &x_face_quad(:,ggi), normal(:,ggi))
+
+                           end if
+
+                           v_dot_n = &
+- di%cached_face_value%relperm(1,ggi,vele,p) * absperm_ele(1) * &
+dot_product((grad_pressure_face_quad(:,ggi) - di%cached_face_value%den(ggi,vele,p) * grav_ele(:,1)), normgi) / &
+visc_ele(1)
+
+                           inflow = (v_dot_n<=0.0)
+
+                           income = merge(1.0,0.0,inflow)
+
+                           cfl_rhs_local(iloc) = cfl_rhs_local(iloc) + &
+                                                 abs(v_dot_n) * &
+                                                 detwei(ggi) * &
+                                                 (1.0 - income)
+
+                           cfl_rhs_local(oloc) = cfl_rhs_local(oloc) + &
+                                                 abs(v_dot_n) * &
+                                                 detwei(ggi) * &
+                                                 income
+
+                           do dim = 1, di%ndim
+
+                              v_local(dim,iloc) = v_local(dim,iloc) - &
+detwei(ggi) * &
+di%cached_face_value%relperm(1,ggi,vele,p) * &
+absperm_ele(1) * &
+(grad_pressure_face_quad(dim,ggi) - di%cached_face_value%den(ggi,vele,p) * grav_ele(dim,1)) / &
+visc_ele(1)
+
+                              v_local(dim,oloc) = v_local(dim,oloc) - &
+detwei(ggi) * &
+di%cached_face_value%relperm(1,ggi,vele,p) * &
+absperm_ele(1) * &
+(grad_pressure_face_quad(dim,ggi) - di%cached_face_value%den(ggi,vele,p) * grav_ele(dim,1)) / &
+visc_ele(1)
+                           
+                           end do
+                           
+                           m_local(iloc) = m_local(iloc) + &
+                                           detwei(ggi) * &
+                                           di%cached_face_value%relperm(1,ggi,vele,p) / &
+                                           visc_ele(1)
+
+                           m_local(oloc) = m_local(oloc) + &
+                                           detwei(ggi) * &
+                                           di%cached_face_value%relperm(1,ggi,vele,p) / &
+                                           visc_ele(1)
+
+                        end if check_visited
+
+                     end do quadrature_loop
+
+                  end if is_neigh
+
+               end do face_loop
+
+            end do nodal_loop_i
+            
+            call addto(di%cfl(p)%ptr, p_nodes, cfl_rhs_local)
+            
+            call addto(di%darcy_velocity(p)%ptr, p_nodes, v_local)
+            
+            call addto(di%mobility(p)%ptr, p_nodes, m_local)
+
+         end do vele_loop
+
+         ! Get this phase v BC info - only for no_normal_flow and
+         ! prescribed_normal_flow which is special as it is a scalar
+         call darcy_impes_get_v_boundary_condition(&
+              di%darcy_velocity(p)%ptr, &
+              (/"prescribed_normal_flow", "no_normal_flow        "/), &
+              di%bc_surface_mesh, di%v_bc_value, di%v_bc_flag)
+         
+         sele_loop: do sele = 1, di%number_sele
+            
+            p_nodes_bdy = face_global_nodes(di%pressure_mesh, sele)
+
+            ! obtain the transformed determinant*weight and normals
+            if (di%cached_face_value%cached_detwei_normal) then
+
+               detwei_bdy => di%cached_face_value%detwei_bdy(:,sele)
+               
+               normal_bdy => di%cached_face_value%normal_bdy(:,:,sele)
+            
+            else
+
+               x_ele     = ele_val(di%positions, face_ele(di%positions, sele))
+               x_ele_bdy = face_val(di%positions, sele)
+            
+               call transform_cvsurf_facet_to_physical(x_ele, x_ele_bdy, di%x_cvbdyshape, normal_bdy, detwei_bdy)
+            
+            end if
+            
+            if (di%v_bc_flag(sele) == V_BC_TYPE_PRESCRIBED_NORMAL_FLOW) then
+            
+               bc_sele_val = ele_val(di%v_bc_value, sele)
+            
+            end if
+            
+            visc_ele_bdy                = face_val(di%viscosity(p)%ptr, sele)
+            absperm_ele_bdy             = face_val(di%absolute_permeability, sele)
+            grad_pressure_face_quad_bdy = face_val_at_quad(di%gradient_pressure(p)%ptr, sele, di%gradp_cvbdyshape)         
+            grav_ele_bdy                = face_val(di%gravity, sele)                         
+            
+            ! Initialise the local assemble arrays for this element
+            cfl_rhs_local_bdy = 0.0
+            v_local_bdy       = 0.0
+            m_local_bdy       = 0.0
+
+            bc_iloc_loop: do iloc = 1, di%pressure_mesh%faces%shape%loc
+
+               bc_face_loop: do face = 1, di%cvfaces%sfaces
+
+                  bc_neigh_if: if(di%cvfaces%sneiloc(iloc,face)/=0) then
+
+                     bc_quad_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                        ggi = (face-1)*di%cvfaces%shape%ngi + gi
+                        
+                        if (di%v_bc_flag(sele) == V_BC_TYPE_NO_NORMAL_FLOW) then
+                        
+                           v_dot_n = 0.0
+                        
+                        else if (di%v_bc_flag(sele) == V_BC_TYPE_PRESCRIBED_NORMAL_FLOW) then
+                        
+                           v_dot_n = bc_sele_val(iloc)
+                                                      
+                        else
+
+                           v_dot_n = &
+- di%cached_face_value%relperm_bdy(1,ggi,sele,p) * absperm_ele_bdy(1) * &
+dot_product((grad_pressure_face_quad_bdy(:,ggi) - di%cached_face_value%den_bdy(ggi,sele,p) * grav_ele_bdy(:,1)), normal_bdy(:,ggi)) / &
+visc_ele_bdy(1)
+                                                                              
+                        end if
+                        
+                        if (v_dot_n > 0.0) then
+                           income = 0.0
+                        else
+                           income = 1.0
+                        end if
+                        
+                        cfl_rhs_local_bdy(iloc) = cfl_rhs_local_bdy(iloc) + &
+                                                  abs(v_dot_n) * &
+                                                  detwei_bdy(ggi) * &
+                                                  (1.0 - income)                     
+
+                        do dim = 1, di%ndim
+
+                           v_local_bdy(dim,iloc) = v_local_bdy(dim,iloc) - &
+detwei_bdy(ggi) * &
+di%cached_face_value%relperm_bdy(1,ggi,sele,p) * &
+absperm_ele_bdy(1) * &
+(grad_pressure_face_quad_bdy(dim,ggi) - di%cached_face_value%den_bdy(ggi,sele,p) * grav_ele_bdy(dim,1)) / &
+visc_ele_bdy(1)
+                           
+                        end do
+
+                        m_local_bdy(iloc) = m_local_bdy(iloc) + &
+                                            detwei_bdy(ggi) * &
+                                            di%cached_face_value%relperm_bdy(1,ggi,sele,p) / &
+                                            visc_ele_bdy(1)
+
+                     end do bc_quad_loop
+
+                  end if bc_neigh_if
+
+               end do bc_face_loop
+
+            end do bc_iloc_loop
+
+            call addto(di%cfl(p)%ptr, p_nodes_bdy, cfl_rhs_local_bdy)
+            
+            call addto(di%darcy_velocity(p)%ptr, p_nodes_bdy, v_local_bdy)
+            
+            call addto(di%mobility(p)%ptr, p_nodes_bdy, m_local_bdy)
+
+         end do sele_loop
+                  
+         di%cfl(p)%ptr%val = di%cfl(p)%ptr%val * di%dt / di%cv_mass_pressure_mesh_with_porosity%val
+         
+         call scale(di%darcy_velocity(p)%ptr, di%inverse_cv_sa_pressure_mesh)
+         
+         call scale(di%mobility(p)%ptr, di%inverse_cv_sa_pressure_mesh)
+         
+         ewrite_minmax(di%cfl(p)%ptr)
+
+         ewrite_minmax(di%darcy_velocity(p)%ptr)
+
+         ewrite_minmax(di%mobility(p)%ptr)
+
+      end do phase_loop      
+            
+      ewrite(1,*) 'Calculate TotalDarcyVelocity and TotalMobility'
+      
+      call set(di%total_darcy_velocity, di%darcy_velocity(1)%ptr)
+      call set(di%total_mobility, di%mobility(1)%ptr)
+      
+      do p = 2, di%number_phase
+         
+         call addto(di%total_darcy_velocity, di%darcy_velocity(p)%ptr)
+         call addto(di%total_mobility, di%mobility(p)%ptr)
+         
+      end do
+
+      ewrite_minmax(di%total_darcy_velocity)
+      ewrite_minmax(di%total_mobility)
+      
+      ewrite(1,*) 'Calculate BulkDarcyVelocity'
+      
+      call allocate(aux_vfield, di%ndim, di%pressure_mesh)
+      
+      call zero(di%bulk_darcy_velocity)
+      
+      do p = 1, di%number_phase
+      
+         call set(aux_vfield, di%darcy_velocity(p)%ptr)
+         call scale(aux_vfield, di%saturation(p)%ptr)
+         
+         call addto(di%bulk_darcy_velocity, aux_vfield)
+         
+      end do
+      
+      call deallocate(aux_vfield)
+      
+      ewrite_minmax(di%bulk_darcy_velocity)
+      
+      ! calculate the fractional flow for each phase
+      do p = 1, di%number_phase
+         
+         ewrite(1,*) 'Calculate FractionalFlow for phase ',p
+                  
+         call set(di%fractional_flow(p)%ptr, di%total_mobility)
+         
+         call invert(di%fractional_flow(p)%ptr)
+         
+         call scale(di%fractional_flow(p)%ptr, di%mobility(p)%ptr)
+                  
+         ewrite_minmax(di%fractional_flow(p)%ptr)
+      
+      end do      
+      
+      ! deallocate local variables as required
+      deallocate(x_ele)
+      if (di%cached_face_value%cached_detwei_normal) then
+         nullify(detwei)
+         nullify(normal)      
+      else
+         deallocate(detwei)
+         deallocate(normal)
+      end if
+      deallocate(normgi)
+      deallocate(notvisited)
+      deallocate(cfl_rhs_local)
+      deallocate(x_face_quad)
+      deallocate(grad_pressure_face_quad)
+      deallocate(grav_ele)
+      
+      deallocate(grav_ele_bdy)
+      deallocate(grad_pressure_face_quad_bdy)
+      if (di%cached_face_value%cached_detwei_normal) then
+         nullify(detwei_bdy)
+         nullify(normal_bdy)      
+      else
+         deallocate(detwei_bdy)
+         deallocate(normal_bdy)
+      end if
+      deallocate(x_ele_bdy)
+      deallocate(p_nodes_bdy)
+      deallocate(cfl_rhs_local_bdy)
+      deallocate(bc_sele_val)      
+
+      deallocate(v_local)
+      deallocate(m_local)
+      deallocate(v_local_bdy)
+      deallocate(m_local_bdy)
+             
+   end subroutine darcy_impes_calculate_vel_mob_ff_and_cfl_fields
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_calculate_inverse_cv_sa(di)
+      
+      !!< Calculate the inverse of the surface area of the CV's on the pressure mesh
+      
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variables
+      integer :: vele, iloc, oloc, face, gi, ggi, sele
+      real,    dimension(:),     pointer :: sa_local
+      real,    dimension(:,:),   pointer :: x_ele
+      real,    dimension(:,:),   pointer :: normal
+      real,    dimension(:),     pointer :: detwei
+      logical, dimension(:),     pointer :: notvisited
+      integer, dimension(:),     pointer :: p_nodes      
+      real,    dimension(:,:),   pointer :: normal_bdy
+      real,    dimension(:),     pointer :: detwei_bdy
+      real,    dimension(:,:),   pointer :: x_ele_bdy
+      real,    dimension(:),     pointer :: sa_local_bdy
+      integer, dimension(:),     pointer :: p_nodes_bdy
+            
+      ewrite(1,*) 'Calculate the inverse of the surface area of CV on the pressure mesh'
+
+      ! allocate arrays used in assemble process - many assume
+      ! that all elements are the same type.
+      allocate(x_ele(di%ndim, ele_loc(di%positions,1)))      
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(normal(di%ndim,di%x_cvshape%ngi))
+         allocate(detwei(di%x_cvshape%ngi))      
+      end if
+      allocate(notvisited(di%x_cvshape%ngi))
+      
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(detwei_bdy(di%x_cvbdyshape%ngi))
+         allocate(normal_bdy(di%ndim, di%x_cvbdyshape%ngi))
+      end if
+      allocate(x_ele_bdy(di%ndim, face_loc(di%positions,1)))      
+      allocate(p_nodes_bdy(face_loc(di%pressure_mesh,1)))      
+      allocate(sa_local(ele_loc(di%pressure_mesh,1)))
+      allocate(sa_local_bdy(face_loc(di%pressure_mesh,1)))
+         
+      call zero(di%inverse_cv_sa_pressure_mesh)
+
+      vele_loop: do vele = 1, di%number_vele
+
+         ! The node indices of the pressure field
+         p_nodes => ele_nodes(di%pressure_mesh, vele)
+
+         ! obtain the transformed determinant*weight and normals
+         if (di%cached_face_value%cached_detwei_normal) then
+
+            detwei => di%cached_face_value%detwei(:,vele)
+
+            normal => di%cached_face_value%normal(:,:,vele)
+
+         else
+
+            ! get the coordinate values for this element for each positions local node
+            x_ele = ele_val(di%positions, vele)         
+
+            call transform_cvsurf_to_physical(x_ele, di%x_cvshape, detwei, normal, di%cvfaces)
+
+         end if
+
+         ! Initialise array for the quadrature points of this 
+         ! element for whether it has already been visited
+         notvisited = .true.
+
+         sa_local = 0.0
+
+         ! loop over local nodes within this element
+         nodal_loop_i: do iloc = 1, di%pressure_mesh%shape%loc
+
+            ! loop over CV faces internal to this element
+            face_loop: do face = 1, di%cvfaces%faces
+
+               ! is this a face neighbouring iloc?
+               is_neigh: if(di%cvfaces%neiloc(iloc, face) /= 0) then
+
+                  ! find the opposing local node across the CV face
+                  oloc = di%cvfaces%neiloc(iloc, face)
+
+                  ! loop over gauss points on face
+                  quadrature_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                     ! global gauss pt index for this element
+                     ggi = (face-1)*di%cvfaces%shape%ngi + gi
+
+                     ! check if this quadrature point has already been visited
+                     check_visited: if(notvisited(ggi)) then
+
+                        notvisited(ggi) = .false.
+
+                        sa_local(iloc) = sa_local(iloc) + &
+                                         detwei(ggi)
+
+                        sa_local(oloc) = sa_local(oloc) + &
+                                         detwei(ggi)
+
+                     end if check_visited
+
+                  end do quadrature_loop
+
+               end if is_neigh
+
+            end do face_loop
+
+         end do nodal_loop_i
+
+         call addto(di%inverse_cv_sa_pressure_mesh, p_nodes, sa_local)
+
+      end do vele_loop
+      
+      sele_loop: do sele = 1, di%number_sele
+
+         p_nodes_bdy = face_global_nodes(di%pressure_mesh, sele)
+
+         ! obtain the transformed determinant*weight and normals
+         if (di%cached_face_value%cached_detwei_normal) then
+
+            detwei_bdy => di%cached_face_value%detwei_bdy(:,sele)
+
+            normal_bdy => di%cached_face_value%normal_bdy(:,:,sele)
+
+         else
+
+            x_ele     = ele_val(di%positions, face_ele(di%positions, sele))
+            x_ele_bdy = face_val(di%positions, sele)
+
+            call transform_cvsurf_facet_to_physical(x_ele, x_ele_bdy, di%x_cvbdyshape, normal_bdy, detwei_bdy)
+
+         end if
+
+         sa_local_bdy = 0.0
+
+         bc_iloc_loop: do iloc = 1, di%pressure_mesh%faces%shape%loc
+
+            bc_face_loop: do face = 1, di%cvfaces%sfaces
+
+               bc_neigh_if: if(di%cvfaces%sneiloc(iloc,face)/=0) then
+
+                  bc_quad_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                     ggi = (face-1)*di%cvfaces%shape%ngi + gi
+
+                     sa_local_bdy(iloc) = sa_local_bdy(iloc) + &
+                                          detwei_bdy(ggi)                    
+
+                  end do bc_quad_loop
+
+               end if bc_neigh_if
+
+            end do bc_face_loop
+
+         end do bc_iloc_loop
+
+         call addto(di%inverse_cv_sa_pressure_mesh, p_nodes_bdy, sa_local_bdy)
+
+      end do sele_loop
+      
+      call invert(di%inverse_cv_sa_pressure_mesh)
+            
+      ! deallocate local variables as required
+      deallocate(x_ele)
+      if (di%cached_face_value%cached_detwei_normal) then
+         nullify(detwei)
+         nullify(normal)      
+      else
+         deallocate(detwei)
+         deallocate(normal)
+      end if
+      deallocate(notvisited)
+      
+      if (di%cached_face_value%cached_detwei_normal) then
+         nullify(detwei_bdy)
+         nullify(normal_bdy)      
+      else
+         deallocate(detwei_bdy)
+         deallocate(normal_bdy)
+      end if
+      deallocate(x_ele_bdy)
+      deallocate(p_nodes_bdy)
+
+      deallocate(sa_local)
+      deallocate(sa_local_bdy)
+      
+      ewrite_minmax(di%inverse_cv_sa_pressure_mesh)
+      
+      ewrite(1,*) 'Finished Calculate the inverse of the surface area of CV on the pressure mesh'
+      
+   end subroutine darcy_impes_calculate_inverse_cv_sa
+
+! ----------------------------------------------------------------------------
+   
+   subroutine darcy_impes_calculate_sum_saturation(di)
+      
+      !!< Calculate the sum of the saturation fields
+      
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variables
+      integer :: p
+      
+      ewrite(1,*) 'Calculate SumSaturation'
+      
+      call zero(di%sum_saturation)
+      
+      do p = 1, di%number_phase
+         
+         call addto(di%sum_saturation, di%saturation(p)%ptr)
+         
+      end do
+      
+      ewrite_minmax(di%sum_saturation)
+       
+   end subroutine darcy_impes_calculate_sum_saturation
+
+! ----------------------------------------------------------------------------
+   
+   subroutine darcy_impes_calculate_densities(di)
+      
+      !!< Calculate the density field of each phase
+      
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variables
+      integer :: p
+      
+      ewrite(1,*) 'Calculate Density of each phase'
+            
+      do p = 1, di%number_phase
+         
+         if (di%eos_options(p)%have_fluids_linear) then
+            
+            ! Currently no variation of the reference density by any pertubation
+            
+            call set(di%density(p)%ptr, di%eos_options(p)%fluids_linear_reference_density)
+         
+         end if
+
+         ewrite_minmax(di%density(p)%ptr)
+         
+      end do 
+
+      ewrite(1,*) 'Finished Calculate Density of each phase'
+       
+   end subroutine darcy_impes_calculate_densities
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_calculate_relperm_fields(di)
+      
+      !!< Calculate the latest relperm fields from the latest saturations.
+      !!< This is assumes that the relperm and saturation fields are 
+      !!< on the same mesh, which is all the input permits
+      
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variables
+      integer                     :: p, node
+      real                        :: relperm_node_val
+      real, dimension(:), pointer :: sat_node_val_all_phases
+
+      ewrite(1,*) 'Calculate RelativePermeability field of each phase'
+      
+      allocate(sat_node_val_all_phases(di%number_phase))
+      
+      ! Calc relperm phase values for each node
+      node_loop: do node = 1, di%number_pmesh_node
+         
+         ! Find all the phase saturation node values
+         do p = 1, di%number_phase
+            
+            sat_node_val_all_phases(p) = node_val(di%saturation(p)%ptr, node)
+            
+         end do
+         
+         ! Calc and set relperm node value for each phase
+         do p = 1, di%number_phase
+            
+            call darcy_impes_calculate_relperm_value(relperm_node_val, &
+                                                     sat_node_val_all_phases, &
+                                                     p, &
+                                                     di%relperm_corr_options%type, &
+                                                     di%relperm_corr_options%exponents, &
+                                                     di%relperm_corr_options%residual_saturations, &
+                                                     di%relperm_corr_options%cutoff_saturations, &
+                                                     di%relperm_corr_options%scaling_coefficients)
+            
+            call set(di%relative_permeability(p)%ptr, &
+                     node, &
+                     relperm_node_val)
+               
+         end do
+         
+      end do node_loop
+
+      deallocate(sat_node_val_all_phases)
+      
+      ewrite(1,*) 'Finished Calculate RelativePermeability field of each phase'
+            
+   end subroutine darcy_impes_calculate_relperm_fields
+
+! ----------------------------------------------------------------------------
+
+   subroutine darcy_impes_calculate_relperm_value(relperm_val, &
+                                                  sat_val_all_phases, &
+                                                  p, &
+                                                  relperm_corr_type, &
+                                                  relperm_corr_exponents, &
+                                                  relperm_corr_residual_sats, &
+                                                  relperm_corr_cutoff_sats, &
+                                                  relperm_corr_scaling_coefficients)
+      
+      !!< Calculate the latest relperm value for phase p for the 
+      !!< given saturation values of all phases using the given options.
+      
+      real,                  intent(out) :: relperm_val
+      real,    dimension(:), intent(in)  :: sat_val_all_phases
+      integer,               intent(in)  :: p
+      integer,               intent(in)  :: relperm_corr_type
+      real,    dimension(:), intent(in)  :: relperm_corr_exponents
+      real,    dimension(:), intent(in)  :: relperm_corr_residual_sats
+      real,    dimension(:), intent(in)  :: relperm_corr_cutoff_sats
+      real,    dimension(:), intent(in)  :: relperm_corr_scaling_coefficients
+    
+      ! local variables
+      real :: sat_minus_res_sat
+      real :: sat_effective
+         
+      select case (relperm_corr_type)
+           
+      case (RELPERM_CORRELATION_POWER)
+
+         relperm_val = (sat_val_all_phases(p) - relperm_corr_residual_sats(p)) ** relperm_corr_exponents(p)
+
+      case (RELPERM_CORRELATION_COREY2PHASE)
+         
+         sat_minus_res_sat = sat_val_all_phases(1) - relperm_corr_residual_sats(1)
+         
+         if (p == 1) then
+
+            relperm_val = sat_minus_res_sat ** 4
+
+         else if (p == 2) then
+
+            relperm_val = (1.0 - sat_minus_res_sat ** 2) * (1.0 - sat_minus_res_sat) ** 2
+
+         else 
+
+            ! This has already been option checked so should not happen
+            FLAbort('Trying to use Corey2Phase relative permeabiltiy correlation for simulation with more than 2 phases')
+
+         end if
+
+      case (RELPERM_CORRELATION_COREY2PHASEOPPOSITE)
+
+         sat_minus_res_sat = sat_val_all_phases(2) - relperm_corr_residual_sats(2)
+
+         if (p == 1) then
+
+            relperm_val = (1.0 - sat_minus_res_sat ** 2) * (1.0 - sat_minus_res_sat) ** 2
+
+         else if (p == 2) then
+
+            relperm_val = sat_minus_res_sat ** 4
+
+         else 
+
+            ! This has already been option checked so should not happen
+            FLAbort('Trying to use Corey2PhaseOpposite relative permeabiltiy correlation for simulation with more than 2 phases')
+
+         end if
+
+      case (RELPERM_CORRELATION_MINERAL)        
+         
+         relperm_val = ((max(relperm_corr_cutoff_sats(p), sat_val_all_phases(p)) - relperm_corr_residual_sats(p)) **&
+                      &relperm_corr_exponents(p)) /  ((1.0 - relperm_corr_residual_sats(p)) ** relperm_corr_exponents(p))
+
+      case (RELPERM_CORRELATION_VANGENUCHTEN)     
+
+         sat_effective = (sat_val_all_phases(2) - relperm_corr_residual_sats(2) ) / ( 1.0 - relperm_corr_residual_sats(1) - &
+                         &relperm_corr_residual_sats(2))
+         
+         if (p == 1) then
+
+            relperm_val = ( 1.0 - sat_effective )**(1.0/3.0) * (1.0 - sat_effective ** ( 1.0/ relperm_corr_exponents(p)) ) ** (2* relperm_corr_exponents(p))
+
+         else if (p == 2) then
+
+            relperm_val =  sat_effective ** (1.0/2.0) * ( 1.0 - ( 1.0 - sat_effective ** ( 1.0 / relperm_corr_exponents(p))) ** relperm_corr_exponents(p))**2.0
+
+         else 
+
+            ! This has already been option checked so should not happen
+            FLAbort('Trying to use VanGenuchten relative permeabiltiy correlation for simulation with more than 2 phases')
+
+         end if         
+
+      case (RELPERM_CORRELATION_JACKSON2PHASE)
+         
+         sat_effective = (sat_val_all_phases(2) - relperm_corr_residual_sats(2)) / &
+                         (1.0 - relperm_corr_residual_sats(1) - relperm_corr_residual_sats(2))
+	 if ( sat_effective >= 1.0 ) then
+	    sat_effective = 1.0
+	 end if 
+         
+	 if ( sat_effective <= 0.0 ) then
+	    sat_effective = 0.0
+	 end if 
+
+         if (p == 1) then
+
+            relperm_val = relperm_corr_scaling_coefficients(1) * (1.0 - sat_effective) ** relperm_corr_exponents(1)
+
+         else if (p == 2) then
+
+            relperm_val = relperm_corr_scaling_coefficients(2) * sat_effective ** relperm_corr_exponents(2)
+
+         else 
+
+            ! This has already been option checked so should not happen
+            FLAbort('Trying to use Jackson2Phase relative permeabiltiy correlation for simulation with more than 2 phases')
+
+         end if
+
+      case (RELPERM_CORRELATION_JACKSON2PHASEOPPOSITE)
+
+         sat_effective = (sat_val_all_phases(1) - relperm_corr_residual_sats(1)) / &
+                         (1.0 - relperm_corr_residual_sats(1) - relperm_corr_residual_sats(2))
+         
+         if (p == 1) then
+
+            relperm_val = relperm_corr_scaling_coefficients(1) * sat_effective ** relperm_corr_exponents(1)
+
+         else if (p == 2) then
+
+            relperm_val = relperm_corr_scaling_coefficients(2) * (1.0 - sat_effective) ** relperm_corr_exponents(2)
+
+         else 
+
+            ! This has already been option checked so should not happen
+            FLAbort('Trying to use Jackson2PhaseOpposite relative permeabiltiy correlation for simulation with more than 2 phases')
+
+         end if
+     
+     end select
+          
+   end subroutine darcy_impes_calculate_relperm_value
+
+! ----------------------------------------------------------------------------
+    
+   subroutine darcy_impes_calculate_cflnumber_field_based_dt(di)
+      
+      !!< Calculate the cfl number field based adaptive time step
+      
+      type(darcy_impes_type), intent(inout) :: di
+
+     !! NEW NoT TESTED FROM HERE
+      type(logical) :: diff_cfl
+
+
+      real :: ds, dpds, diff, deltax
+      type(scalar_field), pointer::sat, phi, perm, vis
+     
+      type(vector_field), pointer:: X
+
+      real, dimension(:), allocatable :: detwei
+      real, dimension(:), pointer :: phi_ele, perm_ele, sat_ele, vis_ele
+      integer, dimension(:), pointer :: sat_ele_nodes
+      integer :: i , ele
+      real, dimension(:, :, :), allocatable :: invJ
+      real :: sat_node
+      type(state_type) :: state
+
+     !! TILL HERE
+  
+      
+      ! local variables
+      integer :: p
+      real, dimension(di%number_phase) :: phase_dt, dtdiff
+      
+      ewrite(1,*) 'Calculate CFL number field based adaptive timestep'
+      
+      ! Find the adaptive dt for each phase then take the minimum
+      
+      phase_dt = di%dt
+      
+      phase_loop: do p = 1, di%number_phase
+        
+         phase_dt(p) = cflnumber_field_based_dt(di%cfl(p)%ptr, &
+                                                di%dt, &
+                                                di%adaptive_dt_options%requested_cfl, &
+                                                di%adaptive_dt_options%min_dt, &
+                                                di%adaptive_dt_options%max_dt, &
+                                                di%adaptive_dt_options%increase_tolerance)         
+         
+      end do phase_loop
+
+
+     ! NEW NOT TESTED FROM HERE 
+!     diff_cfl = .False.        	
+!     if (diff_cfl) then 
+
+!      do p = 1, di%number_phase
+     
+!        sat=>extract_scalar_field(state(p), "Saturation")
+!        phi=>extract_scalar_field(state(p), "Porosity")
+!        perm=>extract_scalar_field(state(p), "Permeability")
+!        vis=>extract_scalar_field(state(p), "Viscosity")
+
+!        X=>extract_vector_field(state(1), "Coordinate")
+!        allocate (invj(mesh_dim(x), mesh_dim(x), ele_ngi(x, 1)))
+!        allocate (detwei(ele_ngi(X,ele)))
+
+
+!        ds = 1.e-10
+!        dpds=0.
+    
+!        diff= 0.
+
+!        do ele = 1, ele_count(sat)
+        
+!          sat_ele_nodes => ele_nodes(sat, ele)
+ 
+!          phi_ele = ele_val(phi, ele)
+!          perm_ele = ele_val(perm, ele)
+!          vis_ele = ele_val(vis, ele)
+
+ !         call compute_inverse_jacobian(ele_val(X,ele), ele_shape(X,ele), &
+  !                                     detwei=detwei, invJ=invJ)
+
+
+ 
+!          deltax = minval(invJ) 
+
+!          do i  = 1, size(sat_ele_nodes)
+ 
+!             sat_node = node_val(sat, sat_ele_nodes(i) )
+!             dpds =  abs( ( (sat_node+ds)**(-0.5) - (sat_node-ds)**(-0.5) ) /(2.*ds) )
+
+ !            diff = max( diff, ( 1./phi_ele(ele)) * (perm_ele(ele) / vis_ele(ele)) * dpds )
+
+!          end do 
+         
+!          dtdiff(p) = deltax**2 / (2*diff)
+
+ !       end do
+
+!       deallocate (invj)
+!       deallocate (detwei)
+
+!     end do
+     
+
+      
+!      di%dt = min(minval(dtdiff), minval(phase_dt))
+
+!     end if 
+
+! UNTIL HERE
+
+      di%dt = minval(phase_dt)
+
+      if(di%adaptive_dt_options%min_dt_terminate_if_reached) then
+         if(di%dt <= di%adaptive_dt_options%min_dt + spacing(di%adaptive_dt_options%min_dt)) then
+            ewrite(0, *) "Minimum timestep reached - terminating"
+            SIG_INT = .true.
+         end if
+      end if
+
+      ewrite(1,*) 'Finished Calculate CFL number field based adaptive timestep'
+       
+   end subroutine darcy_impes_calculate_cflnumber_field_based_dt
+
+! ----------------------------------------------------------------------------
+   
+   subroutine darcy_impes_initialise_cached_face_value(di)
+   
+      !!< Initialise the cached_face_value variable's by 
+      !!< allocating and zeroing the values and setting the 
+      !!< detwei, normal and p_dshape components as required.
+
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variable
+      integer :: vele, sele, iloc, face, gi, ggi
+      real,    dimension(:),   allocatable :: normgi
+      logical, dimension(:),   allocatable :: notvisited
+      real,    dimension(:,:), allocatable :: x_face_quad
+      integer, dimension(:),   pointer     :: x_pmesh_nodes
+
+      ewrite(1,*) 'Initialise cached phase face value types'
+            
+      if (di%subcy_opt_sat%have .and. di%subcy_opt_sat%consistent) then
+      
+         allocate(di%cached_face_value%relperm(di%subcy_opt_sat%number, di%p_cvshape%ngi, di%number_vele, di%number_phase))
+         allocate(di%cached_face_value%relperm_bdy(di%subcy_opt_sat%number, di%p_cvbdyshape%ngi, di%number_sele, di%number_phase))
+         
+      else
+         
+         allocate(di%cached_face_value%relperm(1, di%p_cvshape%ngi, di%number_vele, di%number_phase))
+         allocate(di%cached_face_value%relperm_bdy(1, di%p_cvbdyshape%ngi, di%number_sele, di%number_phase))
+      
+      end if
+      
+      di%cached_face_value%relperm     = 0.0
+      di%cached_face_value%relperm_bdy = 0.0
+      
+      allocate(di%cached_face_value%den(di%p_cvshape%ngi, di%number_vele, di%number_phase))
+      allocate(di%cached_face_value%den_bdy(di%p_cvbdyshape%ngi, di%number_sele, di%number_phase))
+      
+      di%cached_face_value%den     = 0.0
+      di%cached_face_value%den_bdy = 0.0
+      
+      if (di%cached_face_value%cached_detwei_normal) then
+      
+         allocate(di%cached_face_value%detwei(di%p_cvshape%ngi, di%number_vele))
+         allocate(di%cached_face_value%detwei_bdy(di%p_cvbdyshape%ngi, di%number_sele))
+     
+         di%cached_face_value%detwei     = 0.0
+         di%cached_face_value%detwei_bdy = 0.0
+      
+         allocate(di%cached_face_value%normal(di%ndim, di%p_cvshape%ngi, di%number_vele))
+         allocate(di%cached_face_value%normal_bdy(di%ndim, di%p_cvbdyshape%ngi, di%number_sele))
+
+         di%cached_face_value%normal     = 0.0
+         di%cached_face_value%normal_bdy = 0.0
+      
+      end if
+      
+      if (di%cached_face_value%cached_p_dshape) then
+      
+         allocate(di%cached_face_value%p_dshape(di%pressure_mesh%shape%loc, di%p_cvshape%ngi, di%ndim, di%number_vele))
+!!!!! *** THIS IS NOT POSSIBLE YET ***
+!!!!!         allocate(di%cached_face_value%p_dshape_bdy(di%pressure_mesh%faces%shape%loc, di%p_cvbdyshape_full%ngi, di%ndim, di%number_sele))
+
+         di%cached_face_value%p_dshape     = 0.0
+!!!!! *** THIS IS NOT POSSIBLE YET ***
+!!!!!         di%cached_face_value%p_dshape_bdy = 0.0
+      
+      end if
+
+      allocate(normgi(di%ndim))
+      allocate(notvisited(di%p_cvshape%ngi))
+      allocate(x_face_quad(di%ndim, di%p_cvshape%ngi))
+            
+      ! Determine the detwei and normal if required
+      find_cached_detwei_normal: if (di%cached_face_value%cached_detwei_normal) then
+      
+         vele_loop_detwei_normal: do vele = 1, di%number_vele
+
+            call transform_cvsurf_to_physical(ele_val(di%positions, vele), &
+                                              di%x_cvshape, &
+                                              di%cached_face_value%detwei(:,vele), &
+                                              di%cached_face_value%normal(:,:,vele), &
+                                              di%cvfaces)
+            
+            ! Correct the orientation of the normal
+
+            ! get the coordinate values for this element for each quadrature point
+            x_face_quad = ele_val_at_quad(di%positions, vele, di%x_cvshape)
+
+            ! The node indices of the positions projected to the pressure mesh
+            x_pmesh_nodes => ele_nodes(di%positions_pressure_mesh, vele)
+            
+            ! Initialise array for the quadrature points of this 
+            ! element for whether it has already been visited
+            notvisited = .true.
+            
+            ! loop over local nodes within this element
+            nodal_loop_i: do iloc = 1, di%pressure_mesh%shape%loc
+
+              ! loop over CV faces internal to this element
+              face_loop: do face = 1, di%cvfaces%faces
+
+                ! is this a face neighbouring iloc?
+                is_neigh: if(di%cvfaces%neiloc(iloc, face) /= 0) then
+
+                  ! loop over gauss points on face
+                  quadrature_loop: do gi = 1, di%cvfaces%shape%ngi
+
+                    ! global gauss pt index for this element
+                    ggi = (face-1)*di%cvfaces%shape%ngi + gi
+
+                    ! check if this quadrature point has already been visited
+                    check_visited: if(notvisited(ggi)) then
+
+                       notvisited(ggi) = .false.
+
+                       ! correct the orientation of the normal so it points away from iloc
+                       normgi = orientate_cvsurf_normgi(node_val(di%positions_pressure_mesh, x_pmesh_nodes(iloc)), &
+                                                        x_face_quad(:,ggi), &
+                                                        di%cached_face_value%normal(:,ggi,vele))
+                       
+                       di%cached_face_value%normal(:,ggi,vele) = normgi
+                       
+                    end if check_visited
+
+                  end do quadrature_loop
+
+                end if is_neigh
+
+              end do face_loop
+
+            end do nodal_loop_i
+            
+         end do vele_loop_detwei_normal
+
+         sele_loop_detwei_normal: do sele = 1, di%number_sele
+
+            call transform_cvsurf_facet_to_physical(ele_val(di%positions, face_ele(di%positions, sele)), &
+                                                    face_val(di%positions, sele), &
+                                                    di%x_cvbdyshape, &
+                                                    di%cached_face_value%normal_bdy(:,:,sele), &
+                                                    di%cached_face_value%detwei_bdy(:,sele))
+
+         end do sele_loop_detwei_normal
+      
+      end if find_cached_detwei_normal
+
+      deallocate(normgi)
+      deallocate(notvisited)
+      deallocate(x_face_quad)
+      
+      ! Determine the p_dshape if required
+      find_cached_p_dshape: if (di%cached_face_value%cached_p_dshape) then
+
+         vele_loop_p_dshape: do vele = 1, di%number_vele
+            
+            call transform_to_physical(di%positions, &
+                                       vele, &
+                                       x_shape = di%x_cvshape_full, &
+                                       shape   = di%p_cvshape_full, &
+                                       dshape  = di%cached_face_value%p_dshape(:,:,:,vele))
+
+         end do vele_loop_p_dshape
+         
+!!!!! *** THIS IS NOT POSSIBLE YET *** 
+!!!!!         sele_loop_p_dshape: do sele = 1, di%number_sele
+            
+!!!!!            call transform_to_physical(di%positions, &
+!!!!!                                       sele, &
+!!!!!                                       x_shape = di%x_cvbdyshape_full, &
+!!!!!                                       shape   = di%p_cvbdyshape_full, &
+!!!!!                                       dshape  = di%cached_face_value%p_dshape_bdy(:,:,:,sele))
+
+!!!!!         end do sele_loop_p_dshape                 
+         
+      end if find_cached_p_dshape
+      
+      ewrite(1,*) 'Finished Initialise cached face value types'
+      
+   end subroutine darcy_impes_initialise_cached_face_value
+
+! ----------------------------------------------------------------------------
+   
+   subroutine darcy_impes_calculate_relperm_den_first_face_values(di)
+   
+      !!< Find the relperm and density CV first face values for all phase and cache.
+            
+      type(darcy_impes_type), intent(inout) :: di
+      
+      ! local variables
+      logical :: inflow
+      integer :: p, vele, sele, iloc, oloc, face, gi, ggi
+      integer :: r_upwind_pos, d_upwind_pos, s_upwind_pos, dim, sat_p
+      real    :: income, v_over_relperm_dot_n, old_relperm_face_value, den_face_value, old_sat_face_value
+      real,    dimension(1)              :: absperm_ele, visc_ele
+      real,    dimension(:),     pointer :: old_relperm_ele
+      real,    dimension(:),     pointer :: den_ele
+      real,    dimension(:,:),   pointer :: old_sat_ele
+      real,    dimension(:,:),   pointer :: grav_ele
+      real,    dimension(:,:),   pointer :: grad_pressure_face_quad
+      real,    dimension(:,:),   pointer :: v_over_relperm_face_quad
+      real,    dimension(:,:),   pointer :: x_ele
+      real,    dimension(:,:,:), pointer :: p_dshape
+      real,    dimension(:,:),   pointer :: normal
+      real,    dimension(:),     pointer :: detwei
+      real,    dimension(:),     pointer :: normgi
+      logical, dimension(:),     pointer :: notvisited
+      real,    dimension(:,:),   pointer :: x_face_quad
+      integer, dimension(:),     pointer :: x_pmesh_nodes
+      integer, dimension(:),     pointer :: p_nodes      
+      integer, dimension(:),     pointer :: r_upwind_nodes
+      integer, dimension(:),     pointer :: d_upwind_nodes
+      integer, dimension(:),     pointer :: s_upwind_nodes
+      real,    dimension(:),     pointer :: den_ele_bdy
+      real,    dimension(:),     pointer :: old_relperm_ele_bdy
+            
+      ! allocate arrays used in assemble process - many assume
+      ! that all elements are the same type.
+      allocate(x_ele(di%ndim, ele_loc(di%positions,1)))      
+      if (.not. di%cached_face_value%cached_p_dshape) then
+         allocate(p_dshape(ele_loc(di%pressure_mesh,1), di%x_cvshape%ngi, mesh_dim(di%pressure_mesh)))
+      end if
+      if (.not. di%cached_face_value%cached_detwei_normal) then
+         allocate(normal(di%ndim,di%x_cvshape%ngi))
+         allocate(detwei(di%x_cvshape%ngi))      
+      end if
+      allocate(normgi(di%ndim))
+      allocate(notvisited(di%x_cvshape%ngi))
+      allocate(x_face_quad(di%ndim, di%x_cvshape%ngi))
+      allocate(old_relperm_ele(ele_loc(di%pressure_mesh,1)))
+      allocate(den_ele(ele_loc(di%pressure_mesh,1)))
+      allocate(old_sat_ele(ele_loc(di%pressure_mesh,1), di%number_phase))
+      allocate(grav_ele(di%ndim,1))
+      allocate(grad_pressure_face_quad(di%ndim, di%p_cvshape%ngi))
+      allocate(v_over_relperm_face_quad(di%ndim,di%p_cvshape%ngi))
+
+      allocate(den_ele_bdy(face_loc(di%pressure_mesh,1)))
+      allocate(old_relperm_ele_bdy(face_loc(di%pressure_mesh,1)))      
+      
+      ewrite(1,*) 'Calculate relperm and density and first face values'      
+
+      ! Initialise the cached face values
+      di%cached_face_value%relperm(1,:,:,:)     = 0.0
+      di%cached_face_value%relperm_bdy(1,:,:,:) = 0.0
+      di%cached_face_value%den                  = 0.0
+      di%cached_face_value%den_bdy              = 0.0
+
+      phase_loop: do p = 1, di%number_phase
+          
+         ! Determine the upwind relperm values of all node pairs if required for higher order CV face value.
+         ! If the relperm face value is a RELPERMOVERSAT... scheme then the upwind values 
+         ! of modrelperm are required (= relperm / sat)
+         if(di%relperm_cv_options%limit_facevalue) then
+           
+            if (di%determine_saturation_face_values) then
+               
+               call darcy_impes_calculate_old_modified_relative_permeability(di, p)
+               
+               ! NOTE only old values are required but this procedure 
+               !      expects latest and old. So pass in old twice - not optimal
+               call find_upwind_values(di%state, &
+                                       di%positions_pressure_mesh, &
+                                       di%modified_relative_permeability, &
+                                       di%relperm_upwind, &
+                                       di%modified_relative_permeability, &
+                                       di%relperm_upwind, &
+                                       option_path = trim(di%relative_permeability(p)%ptr%option_path))
+           
+            else 
+           
+               ! NOTE only old values are required but this procedure 
+               !      expects latest and old. So pass in old twice - not optimal
+               call find_upwind_values(di%state, &
+                                       di%positions_pressure_mesh, &
+                                       di%old_relative_permeability(p)%ptr, &
+                                       di%relperm_upwind, &
+                                       di%old_relative_permeability(p)%ptr, &
+                                       di%relperm_upwind, &
+                                       option_path = trim(di%relative_permeability(p)%ptr%option_path))
+
+            end if
+            
+         end if
+
+         ! If the relperm face value is a RELPERMOVERSAT... scheme then the saturation 
+         ! face values need determining which may require the upwind values. 
+         if (di%determine_saturation_face_values .and. di%saturation_cv_options%limit_facevalue) then
+
+            ! NOTE only old values are required but this procedure       
+            !      expects latest and old. So pass in old twice - not optimal
+            call find_upwind_values(di%state, &
+                                    di%positions_pressure_mesh, &
+                                    di%old_saturation(p)%ptr, &
+                                    di%sfield_upwind, &
+                                    di%old_saturation(p)%ptr, &
+                                    di%sfield_upwind, &
+                                    option_path = trim(di%saturation(p)%ptr%option_path))            
+
+         end if 
+                          
+         ! Initialise optimisation flag used in finding upwind value in high resolution schemes
+         r_upwind_pos = 0
+         s_upwind_pos = 0
+            
+         vol_element_loop_relperm: do vele = 1, di%number_vele
+           
+            ! The node indices of the pressure mesh
+            p_nodes => ele_nodes(di%pressure_mesh, vele)
+
+            ! The node indices of the positions projected to the pressure mesh
+            x_pmesh_nodes => ele_nodes(di%positions_pressure_mesh, vele)
+                        
+            ! Determine the node numbers to use to determine the relperm upwind values
+            if((di%relperm_cv_options%upwind_scheme == CV_UPWINDVALUE_PROJECT_POINT).or.&
+               (di%relperm_cv_options%upwind_scheme == CV_UPWINDVALUE_PROJECT_GRAD)) then
+
+               r_upwind_nodes => x_pmesh_nodes
+
+            else
+
+               r_upwind_nodes => p_nodes
+
+            end if
+
+            if (di%determine_saturation_face_values) then
+            
+               ! Determine the node numbers to use to determine the saturation upwind values
+               if((di%saturation_cv_options%upwind_scheme == CV_UPWINDVALUE_PROJECT_POINT).or.&
+                  (di%saturation_cv_options%upwind_scheme == CV_UPWINDVALUE_PROJECT_GRAD)) then
+
+                  s_upwind_nodes => x_pmesh_nodes
+
+               else
+
+                  s_upwind_nodes => p_nodes
+
+               end if
+               
+               do sat_p = 1, di%number_phase
+               
+                  ! get the old saturation ele values for this phase
+                  old_sat_ele(:,sat_p) = ele_val(di%old_saturation(sat_p)%ptr, vele)
+               
+               end do
+               
+            end if
+
+            ! get the old_relperm ele values for this phase
+            old_relperm_ele = ele_val(di%old_relative_permeability(p)%ptr, vele)
+            
+            ! The gravity values for this element for each direction
+            grav_ele = ele_val(di%gravity, vele) 
+
+            ! get the viscosity value for this element for this phase
+            visc_ele = ele_val(di%viscosity(p)%ptr, vele)
+
+            ! get the absolute permeability value for this element
+            absperm_ele = ele_val(di%absolute_permeability, vele)         
+            
+            ! obtain the transformed determinant*weight and normals
+            if (di%cached_face_value%cached_detwei_normal) then
+               
+               detwei => di%cached_face_value%detwei(:,vele)
+               
+               normal => di%cached_face_value%normal(:,:,vele)
+               
+            else
+
+               ! get the coordinate values for this element for each positions local node
+               x_ele = ele_val(di%positions, vele)         
+
+               ! get the coordinate values for this element for each quadrature point
+               x_face_quad = ele_val_at_quad(di%positions, vele, di%x_cvshape)
+            
+               call transform_cvsurf_to_physical(x_ele, di%x_cvshape, detwei, normal, di%cvfaces)
+            
+            end if
+            
+            ! obtain the derivative of the pressure mesh shape function at the CV face quadrature points
+            if (di%cached_face_value%cached_p_dshape) then
+               
+               p_dshape => di%cached_face_value%p_dshape(:,:,:,vele)
+               
+            else
+            
+               call transform_to_physical(di%positions, vele, x_shape = di%x_cvshape_full, &
+                                          shape = di%p_cvshape_full, dshape = p_dshape)
+            
+            end if
+            
+            ! get the gradient pressure at the cv surface quadrature points for each direction for this phase
+            grad_pressure_face_quad = darcy_impes_ele_grad_at_quad_scalar(di%pressure(p)%ptr, vele, dn = p_dshape)
+                        
+            ! the latest DarcyVelocity/relperm at the quadrature points
+            ! determined from FE interpolation of each component, only used to determine upwind.
+            do dim = 1,di%ndim
+
+               v_over_relperm_face_quad(dim,:) = - (absperm_ele(1) / visc_ele(1)) * &
+(grad_pressure_face_quad(dim,:) - ele_val_at_quad(di%density(p)%ptr, vele, di%p_cvshape) * grav_ele(dim,1))
+
+            end do
+            
+            ! Initialise array for the quadrature points of this 
+            ! element for whether it has already been visited
+            notvisited = .true.
+            
+            ! loop over local nodes within this element
+            nodal_loop_i_relperm: do iloc = 1, di%pressure_mesh%shape%loc
+
+              ! loop over CV faces internal to this element
+              face_loop_relperm: do face = 1, di%cvfaces%faces
+
+                ! is this a face neighbouring iloc?
+                is_neigh_relperm: if(di%cvfaces%neiloc(iloc, face) /= 0) then
+
+                  ! find the opposing local node across the CV face
+                  oloc = di%cvfaces%neiloc(iloc, face)
+
+                  ! loop over gauss points on face
+                  quadrature_loop_relperm: do gi = 1, di%cvfaces%shape%ngi
+
+                    ! global gauss pt index for this element
+                    ggi = (face-1)*di%cvfaces%shape%ngi + gi
+
+                    ! check if this quadrature point has already been visited
+                    check_visited_relperm: if(notvisited(ggi)) then
+
+                       notvisited(ggi) = .false.
+
+                       if (di%cached_face_value%cached_detwei_normal) then
+                          
+                          ! cached normal has correct orientation already
+                          normgi = normal(:,ggi)
+                          
+                       else
+                       
+                          ! correct the orientation of the normal so it points away from iloc
+                          normgi = orientate_cvsurf_normgi(node_val(di%positions_pressure_mesh, x_pmesh_nodes(iloc)), &
+                                                          &x_face_quad(:,ggi), normal(:,ggi))
+                       
+                       end if
+                       
+                       ! determine if the flow is in or out of the face at this quadrature
+                       ! with respect to the normal orientation using the latest v_over_relperm
+                       v_over_relperm_dot_n = dot_product(v_over_relperm_face_quad(:,ggi), normgi)
+
+                       inflow = (v_over_relperm_dot_n<=0.0)
+
+                       income = merge(1.0,0.0,inflow)
+
+                       ! Evaluate the face value for relperm 
+                       call darcy_impes_evaluate_relperm_face_val(old_relperm_face_value, &
+                                                                  iloc, &
+                                                                  oloc, &
+                                                                  ggi, &
+                                                                  r_upwind_nodes, &
+                                                                  di%p_cvshape, &
+                                                                  old_relperm_ele, &
+                                                                  di%relperm_upwind, &
+                                                                  inflow, &
+                                                                  income, &
+                                                                  old_sat_ele, &
+                                                                  di%minimum_denominator_saturation_value, &
+                                                                  p, &
+                                                                  di%number_phase, &
+                                                                  di%relperm_corr_options, &
+                                                                  di%relperm_cv_options, &
+                                                                  save_pos = r_upwind_pos)
+
+                       ! Evaluate the face value for saturation if required, else set to 1.0
+                       if (di%determine_saturation_face_values) then
+                       
+                          call darcy_impes_evaluate_sfield_face_val(old_sat_face_value, &
+                                                                    iloc, &
+                                                                    oloc, &
+                                                                    ggi, &
+                                                                    s_upwind_nodes, &
+                                                                    di%p_cvshape, &
+                                                                    old_sat_ele(:,p), &
+                                                                    di%sfield_upwind, &
+                                                                    inflow, &
+                                                                    income, &
+                                                                    di%saturation_cv_options, &
+                                                                    save_pos = s_upwind_pos)                          
+                       
+                       else
+                       
+                          old_sat_face_value = 1.0
+                                              
+                       end if
+                                              
+                       ! Cache the old relperm (or modrelperm*sat) face value
+                       di%cached_face_value%relperm(1,ggi,vele,p) = old_relperm_face_value * old_sat_face_value
+                                           
+                    end if check_visited_relperm
+
+                  end do quadrature_loop_relperm
+
+                end if is_neigh_relperm
+
+              end do face_loop_relperm
+
+            end do nodal_loop_i_relperm
+
+         end do vol_element_loop_relperm
+ 
+         ! Determine the upwind density values of all node pairs if required for higher order CV face value
+         if(di%density_cv_options%limit_facevalue) then
+
+           ! NOTE only new upwind values are required but this procedure 
+           !      expects latest and old. So pass in new twice - not optimal
+           call find_upwind_values(di%state, &
+                                   di%positions_pressure_mesh, &
+                                   di%density(p)%ptr, &
+                                   di%sfield_upwind, &
+                                   di%density(p)%ptr, &
+                                   di%sfield_upwind, &
+                                   option_path = trim(di%density(p)%ptr%option_path))
+
+         end if
+                          
+         ! Initialise optimisation flag used in finding upwind value in high resolution schemes
+         d_upwind_pos = 0
+            
+         vol_element_loop_den: do vele = 1, di%number_vele
+           
+            ! The node indices of the pressure mesh
+            p_nodes => ele_nodes(di%pressure_mesh, vele)
+
+            ! The node indices of the positions projected to the pressure mesh
+            x_pmesh_nodes => ele_nodes(di%positions_pressure_mesh, vele)
+            
+            ! Determine the node numbers to use to determine the density upwind values
+            if((di%density_cv_options%upwind_scheme == CV_UPWINDVALUE_PROJECT_POINT).or.&
+               (di%density_cv_options%upwind_scheme == CV_UPWINDVALUE_PROJECT_GRAD)) then
+
+               d_upwind_nodes => x_pmesh_nodes
+
+            else
+
+               d_upwind_nodes => p_nodes
+
+            end if
+
+            ! get the density value for this element for this phase
+            den_ele = ele_val(di%density(p)%ptr, vele)
+            
+            ! The gravity values for this element for each direction
+            grav_ele = ele_val(di%gravity, vele) 
+
+            ! get the viscosity value for this element for this phase
+            visc_ele = ele_val(di%viscosity(p)%ptr, vele)
+
+            ! get the absolute permeability value for this element
+            absperm_ele = ele_val(di%absolute_permeability, vele)         
+            
+            ! obtain the transformed determinant*weight and normals
+            if (di%cached_face_value%cached_detwei_normal) then
+               
+               detwei => di%cached_face_value%detwei(:,vele)
+               
+               normal => di%cached_face_value%normal(:,:,vele)
+               
+            else
+
+               ! get the coordinate values for this element for each positions local node
+               x_ele = ele_val(di%positions, vele)         
+
+               ! get the coordinate values for this element for each quadrature point
+               x_face_quad = ele_val_at_quad(di%positions, vele, di%x_cvshape)
+            
+               call transform_cvsurf_to_physical(x_ele, di%x_cvshape, detwei, normal, di%cvfaces)
+            
+            end if
+            
+            ! obtain the derivative of the pressure mesh shape function at the CV face quadrature points
+            if (di%cached_face_value%cached_p_dshape) then
+               
+               p_dshape => di%cached_face_value%p_dshape(:,:,:,vele)