~reducedmodelling/fluidity/ROM_Non-intrusive-ann

!    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

!

!    amcgsoftware@imperial.ac.uk

!    

!    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 field_equations_cv

  !!< This module contains the assembly subroutines for advection

  !!< using control volumes

  use fields

  use sparse_matrices_fields

  use sparsity_patterns_meshes

  use state_module

  use spud

  use cv_shape_functions

  use cv_faces

  use cvtools

  use cv_fields

  use cv_upwind_values

  use cv_face_values

  use diagnostic_fields, only: calculate_diagnostic_variable

  use cv_options

  use diagnostic_variables, only: field_tag

  use boundary_conditions

  use boundary_conditions_from_options

  use divergence_matrix_cv, only: assemble_divergence_matrix_cv

  use global_parameters, only: OPTION_PATH_LEN

  use fefields, only: compute_cv_mass

  use petsc_solve_state_module

  use transform_elements, only: transform_cvsurf_to_physical, &

                                transform_cvsurf_facet_to_physical

  use parallel_tools, only: getprocno

  use halos

  use field_options

  use state_fields_module

  use porous_media

  implicit none

  private

  public :: solve_field_eqn_cv, field_equations_cv_check_options, &

            initialise_advection_convergence, coupled_cv_field_eqn, &

            assemble_advectiondiffusion_m_cv

  integer, dimension(:), allocatable, save :: conv_unit

  !! This allows a reference between the field and the file its meant to be

  !! writing to.

  character(len=OPTION_PATH_LEN), dimension(:), allocatable, save :: sfield_list

  ! are we moving the mesh?

  logical :: move_mesh = .false.

  ! are we including density?

  logical :: include_density = .false.

  ! are we including a souce?

  logical :: include_source = .false.

  ! Add source directly to the right hand side?

  logical :: add_src_directly_to_rhs = .false.

  ! are we including an absorption?

  logical :: include_absorption = .false.

  ! are we including diffusion?

  logical :: include_diffusion = .false.

  ! are we including advection?

  logical :: include_advection = .true.

  ! are we including mass?

  logical :: include_mass = .true.

  ! are we assembling particular matrices?

  ! advection?

  logical :: assemble_advection_matrix = .true.

  ! diffusion?

  logical :: assemble_diffusion = .false.

contains

    !************************************************************************

    ! solution wrapping subroutines

    subroutine solve_field_eqn_cv(field_name, state, global_it)

      !!< Construct and solve the advection-diffusion equation for the given

      !!< field using control volumes.

      !! Name of the field to be solved for.

      character(len=*), intent(in) :: field_name

      !! Collection of fields defining system state.

      type(state_type), dimension(:), intent(inout) :: state

      ! global iteration number - passed in so it can be output to file

      integer, intent(in) :: global_it

      ! Field to be solved for (plus its old and globally iterated versions)

      type(scalar_field), pointer :: tfield, oldtfield, it_tfield

      ! Density fields associated with tfield's equation (i.e. if its not pure advection)

      type(scalar_field), pointer :: tdensity, oldtdensity

      ! Coordinate field(s)

      type(vector_field), pointer :: x, x_old, x_new

      type(vector_field) :: x_tfield

      type(tensor_field), pointer :: diffusivity

      type(scalar_field), pointer :: source, absorption

      ! Change in tfield over one timestep.

      type(scalar_field) :: delta_tfield

      ! LHS equation matrix.

      type(csr_matrix) :: M

      ! Advection matrix

      type(csr_matrix) :: A_m

      ! Diffusion matrix

      type(csr_matrix) :: D_m

      ! sparsity structure to construct the matrices with

      type(csr_sparsity), pointer :: mesh_sparsity_1, mesh_sparsity, &

                                     mesh_sparsity_x, grad_m_t_sparsity

      ! Right hand side vector, cv mass matrix, 

      ! locally iterated field (for advection iterations) 

      ! and local old field (for subcycling)

      type(scalar_field), pointer :: t_cvmass, q_cvmass, t_abs_src_cvmass

      type(scalar_field) :: t_cvmass_old, t_cvmass_new

      type(scalar_field) :: rhs, cvmass, advit_tfield, l_old_tfield

      ! Diffusion contribution to rhs

      type(scalar_field) :: diff_rhs

      ! Porosity field

      type(scalar_field) :: porosity_theta 

      type(scalar_field), target :: t_cvmass_with_porosity

      logical :: include_porosity

      ! local copy of option_path for solution field

      character(len=OPTION_PATH_LEN) :: option_path, tdensity_option_path

      ! number of advection iterations and subcycles

      integer :: adv_iterations, no_subcycles

      ! iterators

      integer :: adv_it, sub

      ! time (to output to file), timestep, iterations tolerance, subcycling timestep

      real :: time, dt, error, adv_tolerance, sub_dt

      ! degree of quadrature to use on each control volume face

      integer :: quaddegree

      ! control volume face information

      type(cv_faces_type) :: cvfaces

      ! control volume shape function for volume and boundary

      type(element_type) :: u_cvshape, u_cvbdyshape

      type(element_type) :: ug_cvshape, ug_cvbdyshape

      type(element_type) :: x_cvshape, x_cvbdyshape, &

                            x_cvshape_full, x_cvbdyshape_full

      ! t_cvshape is the element with reduced numbers of derivatives

      ! taken across the control volume faces

      ! t_cvshape_full contains the derivatives with respect to the parent

      ! elements canonical coordinates evaluated at the control volume faces

      ! t_cvbdyshape_full is the same but on the boundary

      type(element_type) :: t_cvshape, t_cvshape_full, diff_cvshape_full, &

                                       t_cvbdyshape_full, diff_cvbdyshape_full

      ! options wrappers for tfield and tdensity

      type(cv_options_type) :: tfield_options, tdensity_options

      ! a dummy density in case we're solving for Advection

      type(scalar_field), pointer :: dummyscalar

      type(vector_field), pointer :: dummyvector

      type(tensor_field), pointer :: dummytensor

      ! somewhere to put strings temporarily

      character(len=FIELD_NAME_LEN) :: tmpstring

      ! what equation type are we solving for?

      integer :: equation_type

      ! success indicator

      integer :: stat

      ! the courant number field

      type(scalar_field) :: cfl_no

      ! nonlinear and grid velocities

      type(vector_field), pointer :: nu, ug

      ! advection velocity

      type(vector_field), pointer :: advu

      !! Gravitational sinking term

      type(scalar_field), pointer :: sink

      !! Direction of gravity

      type(vector_field), pointer :: gravity

      ! assume explicitness?

      logical :: explicit

      ! if we're subcycling how fast can we go?

      real :: max_sub_cfl, max_cfl

      ! temporary hack to get around compiler failure to construct arrays of characters

      character(len=OPTION_PATH_LEN), dimension(1) :: option_path_array

      character(len=OPTION_PATH_LEN), dimension(1) :: tdensity_option_path_array

      ewrite(2,*) 'in solve_field_eqn_cv'

      ewrite(2,*) 'solving for '//field_name//' in state '//trim(state(1)%name)

      ! extract lots of fields:

      ! the actual thing we're trying to solve for

      tfield=>extract_scalar_field(state(1), trim(field_name))

      ewrite_minmax(tfield)

      option_path=tfield%option_path

      ! its previous timelevel - if we're doing more than one iteration per timestep!

      oldtfield=>extract_scalar_field(state(1), "Old"//trim(field_name))

      ! because fluidity resets tfield to oldtfield at the start of every

      ! global iteration we need to undo this so that the control volume faces

      ! are discretised using the most up to date values

      ! therefore extract the iterated values:

      it_tfield=>extract_scalar_field(state(1), "Iterated"//trim(field_name))

      ! and set tfield to them:

      call set(tfield, it_tfield)

      ewrite_minmax(tfield)

      ewrite_minmax(oldtfield)

      ! allocate dummy scalar in case density/source/absorption fields aren't needed (this can be a constant field!)

      allocate(dummyscalar)

      call allocate(dummyscalar, tfield%mesh, name="DummyScalar", field_type=FIELD_TYPE_CONSTANT)

      call zero(dummyscalar)

      dummyscalar%option_path = " "

      ! allocate dummy vector in case the velocity field isn't needed (this can be a constant field!)

      allocate(dummyvector)

      call allocate(dummyvector, mesh_dim(tfield), tfield%mesh, name="DummyVector", field_type=FIELD_TYPE_CONSTANT)

      call zero(dummyvector)

      dummyvector%option_path = " "

      ! allocate dummy tensor in case diffusivity field isn't needed (this can be a constant field!)

      allocate(dummytensor)

      call allocate(dummytensor, tfield%mesh, name="DummyTensor", field_type=FIELD_TYPE_CONSTANT)

      call zero(dummytensor)

      dummytensor%option_path = " "

      ! find out equation type and hence if density is needed or not

      equation_type=equation_type_index(trim(option_path))

      include_density = .false.

      select case(equation_type)

      case(FIELD_EQUATION_ADVECTIONDIFFUSION)

        ! density not needed so use a constant field for assembly

        tdensity=>dummyscalar

        oldtdensity=>dummyscalar

      case(FIELD_EQUATION_CONSERVATIONOFMASS, FIELD_EQUATION_REDUCEDCONSERVATIONOFMASS, &

           FIELD_EQUATION_INTERNALENERGY, FIELD_EQUATION_HEATTRANSFER )

        call get_option(trim(option_path)//'/prognostic/equation[0]/density[0]/name', &

                        tmpstring)

        include_density = .true.

        ! density needed so extract the type specified in the input

        ! ?? are there circumstances where this should be "Iterated"... need to be

        ! careful with priority ordering

        tdensity=>extract_scalar_field(state(1), trim(tmpstring))

        ewrite_minmax(tdensity)

        ! halo exchange? - not currently necessary when suboptimal halo exchange if density

        ! is solved for with this subroutine and the correct priority ordering.

        oldtdensity=>extract_scalar_field(state(1), "Old"//trim(tmpstring))

        ewrite_minmax(oldtdensity)

      end select

      ! get the density option path

      if(have_option(trim(option_path)//'/prognostic/equation[0]/density[0]/discretisation_options')) then

        tdensity_option_path = trim(option_path)//'/prognostic/equation[0]/density[0]/discretisation_options'

      else

        tdensity_option_path = trim(tdensity%option_path)

      end if

      ! now we can get the options for these fields

      ! handily wrapped in a new type...

      tfield_options=get_cv_options(tfield%option_path, tfield%mesh%shape%numbering%family, mesh_dim(tfield))

      if(include_density) then

        tdensity_options=get_cv_options(tdensity_option_path, tdensity%mesh%shape%numbering%family, mesh_dim(tdensity),  coefficient_field=.true.)

      end if

      ! extract fields from state

      ! the base Coordinate field

      x=>extract_vector_field(state(1), "Coordinate")

      ! the Coordinate field on the same mesh as tfield

      x_tfield=get_coordinate_field(state(1), tfield%mesh)

      ! are we including the mass (generally yes)?

      include_mass = .not.have_option(trim(tfield%option_path)//"/prognostic/spatial_discretisation/control_volumes/mass_terms/exclude_mass_terms")

      ! are we inclulding advection (generally yes)?

      include_advection = .not.(tfield_options%facevalue==CV_FACEVALUE_NONE)

      if(include_advection) then

        nu=>extract_vector_field(state(1), "NonlinearVelocity")

        ewrite_minmax(nu)

        ! find relative velocity

        allocate(advu)

        call allocate(advu, nu%dim, nu%mesh, "AdvectionVelocity")

        call set(advu, nu)

        if (have_option(trim(tfield%option_path)// & 

         "/prognostic/spatial_discretisation/control_volumes/"// &

         "face_value::FiniteElement/only_sinking_velocity")) then

            call zero(advu)

            ewrite(2,*) "Removing velocity from ", trim(field_name)

        end if

        ! add in sinking velocity

        sink=>extract_scalar_field(state(1), trim(field_name)//"SinkingVelocity"&

            &, stat=stat)

        if(stat==0) then

          gravity=>extract_vector_field(state(1), "GravityDirection")

          ! this may perform a "remap" internally from CoordinateMesh to VelocityMesh

          call addto(advu, gravity, scale=sink)

        end if

        ewrite_minmax(advu)

      else

        ewrite(2,*) 'Excluding advection'

        advu => dummyvector

        if(has_scalar_field(state(1), trim(field_name)//"SinkingVelocity")) then

            ewrite(-1,*) "No advection in "//trim(field_name)

            FLExit("But you have a sinking velocity on. Can't have that")

        end if

      end if

      ! do we have a diffusivity - this will control whether we construct an auxiliary

      ! eqn or not (if BassiRebay is selected) or whether we construct gradients (if ElementGradient

      ! is selected)

      diffusivity=>extract_tensor_field(state(1), trim(field_name)//"Diffusivity", stat=stat)

      include_diffusion = (stat==0).and.(tfield_options%diffusionscheme/=CV_DIFFUSION_NONE)

      if(.not.include_diffusion) then

        diffusivity => dummytensor

      else

        ewrite_minmax(diffusivity)

      end if

      ! do we have a source?

      source=>extract_scalar_field(state(1), trim(field_name)//"Source", stat=stat)

      include_source = (stat==0)

      if(.not.include_source) then

        source=>dummyscalar

      else

        add_src_directly_to_rhs = have_option(trim(source%option_path)//'/diagnostic/add_directly_to_rhs')

        if (add_src_directly_to_rhs) then 

          ewrite(2, *) "Adding Source field directly to the right hand side"

          assert(node_count(source) == node_count(tfield))

        end if

        ewrite_minmax(source)

      end if

      ! do we have an absorption?

      absorption=>extract_scalar_field(state(1), trim(field_name)//"Absorption", stat=stat)

      include_absorption = (stat==0)

      if(.not.include_absorption) then

        absorption=>dummyscalar

      else

        ewrite_minmax(absorption)

      end if

      ! create control volume shape functions

      call get_option("/geometry/quadrature/controlvolume_surface_degree", &

                     quaddegree, default=1)

      cvfaces=find_cv_faces(vertices=ele_vertices(tfield, 1), &

                            dimension=mesh_dim(tfield), &

                            polydegree=tfield%mesh%shape%degree, &

                            quaddegree=quaddegree)

      x_cvshape=make_cv_element_shape(cvfaces, x%mesh%shape)

      t_cvshape=make_cv_element_shape(cvfaces, tfield%mesh%shape)

      if(include_advection) then

        u_cvshape=make_cv_element_shape(cvfaces, nu%mesh%shape)

      else

        u_cvshape=t_cvshape

        call incref(u_cvshape)

      end if

      x_cvbdyshape=make_cvbdy_element_shape(cvfaces, x%mesh%faces%shape)

      if(include_advection) then

        u_cvbdyshape=make_cvbdy_element_shape(cvfaces, nu%mesh%faces%shape)

      else

        u_cvbdyshape=x_cvbdyshape

        call incref(u_cvbdyshape)

      end if

      if(include_diffusion.and.(tfield_options%diffusionscheme==CV_DIFFUSION_ELEMENTGRADIENT)) then

        x_cvshape_full=make_cv_element_shape(cvfaces, x%mesh%shape, &

                                        type=ELEMENT_CONTROLVOLUME_SURFACE_BODYDERIVATIVES)

        t_cvshape_full=make_cv_element_shape(cvfaces, tfield%mesh%shape, &

                                        type=ELEMENT_CONTROLVOLUME_SURFACE_BODYDERIVATIVES)

        diff_cvshape_full=make_cv_element_shape(cvfaces, diffusivity%mesh%shape, &

                                        type=ELEMENT_CONTROLVOLUME_SURFACE_BODYDERIVATIVES)

        x_cvbdyshape_full=make_cvbdy_element_shape(cvfaces, x%mesh%faces%shape, &

                                        type=ELEMENT_CONTROLVOLUME_SURFACE_BODYDERIVATIVES)

        t_cvbdyshape_full=make_cvbdy_element_shape(cvfaces, tfield%mesh%shape, &

                                        type=ELEMENT_CONTROLVOLUME_SURFACE_BODYDERIVATIVES)

        diff_cvbdyshape_full=make_cvbdy_element_shape(cvfaces, diffusivity%mesh%shape, &

                                        type=ELEMENT_CONTROLVOLUME_SURFACE_BODYDERIVATIVES)

      else

        x_cvshape_full=x_cvshape

        t_cvshape_full=t_cvshape

        diff_cvshape_full=t_cvshape

        x_cvbdyshape_full=x_cvbdyshape

        t_cvbdyshape_full=x_cvbdyshape

        diff_cvbdyshape_full=x_cvbdyshape

        call incref(x_cvshape_full)

        call incref(t_cvshape_full)

        call incref(diff_cvshape_full)

        call incref(x_cvbdyshape_full)

        call incref(t_cvbdyshape_full)

        call incref(diff_cvbdyshape_full)

      end if

      ! is this explicit?

      explicit=have_option(trim(option_path)//"/prognostic/explicit")

      ! find the timestep

      call get_option("/timestepping/timestep", dt)

      call get_option("/timestepping/current_time", time) ! so it can be output in the convergence file

      ! allocate and retrieve the cfl no. if necessary

      option_path_array(1) = trim(option_path)                  ! temporary hack for compiler failure

      tdensity_option_path_array(1) = tdensity_option_path

      call cv_disc_get_cfl_no(option_path_array, &

                      state(1), tfield%mesh, cfl_no, &

                      tdensity_option_path_array)

      ! get the mesh sparsity for the matrices

      if(include_diffusion.and.(tfield_options%diffusionscheme==CV_DIFFUSION_BASSIREBAY)) then

        ! in this case we need to extend the sparsity to second order

        mesh_sparsity => get_csr_sparsity_secondorder(state, tfield%mesh, diffusivity%mesh)

        ! except for some things we still need a first order sparsity (upwind value matrices for instance)

        ! although note that this may get modified when periodic depending on the face value scheme

        mesh_sparsity_1 => get_csr_sparsity_firstorder(state, tfield%mesh, tfield%mesh)

        if(.not.(tfield%mesh==diffusivity%mesh)) then

          if(tfield%mesh%shape%degree>1) then

            FLExit("To have a different diffusivity mesh the field must be at most P1")

          elseif(diffusivity%mesh%shape%degree>tfield%mesh%shape%degree) then

            FLExit("The diffusivity mesh must be of a lower degree than the field")

          end if

          ! we also need a first order sparsity using the diffusivity mesh for the assembly of the auxilliary eqn

          grad_m_t_sparsity => get_csr_sparsity_firstorder(state, tfield%mesh, diffusivity%mesh)

        else

          ! no difference between the meshes so this is safe (and faster)

          grad_m_t_sparsity => mesh_sparsity_1

        end if

      else

        ! no BassiRebay diffusion so only a first order sparsity is needed... woo!

        mesh_sparsity => get_csr_sparsity_firstorder(state, tfield%mesh, tfield%mesh)

        grad_m_t_sparsity => mesh_sparsity

        mesh_sparsity_1 => mesh_sparsity

      end if

      if(mesh_periodic(tfield)) then

        ! we have a periodic mesh and depending on the upwind value scheme

        ! we may want to modify the sparsity for the upwind value matrices

        if((tfield_options%upwind_scheme==CV_UPWINDVALUE_PROJECT_POINT).or.&

           (tfield_options%upwind_scheme==CV_UPWINDVALUE_PROJECT_GRAD)) then

           ! yup, we need an unperiodic sparsity

           mesh_sparsity_x => get_csr_sparsity_firstorder(state, x_tfield%mesh, x_tfield%mesh)

        else

           ! periodic sparsity is fine

           mesh_sparsity_x => mesh_sparsity_1

        end if

      else

        ! not periodic so who cares

        mesh_sparsity_x => mesh_sparsity_1

      end if

      if(.not.explicit) then

        ! allocate the lhs matrix

        call allocate(M, mesh_sparsity, name=trim(field_name)//"Matrix")

        call zero(M)

        ! allocate the advection matrix

        call allocate(A_m, mesh_sparsity, name=trim(field_name)//"AdvectionMatrix")

        call zero(A_m)

      else

        if(.not.include_mass) then

          FLExit("Can't be explicit and exclude the mass terms.")

        end if

        ! allocate a local cvmass field because it will get modified by bc etc.

        call allocate(cvmass, tfield%mesh, name=trim(field_name)//"CVMass")

        call zero(cvmass)

      end if

      if(include_diffusion) then

        call allocate(D_m, sparsity=mesh_sparsity, name=trim(field_name)//"AuxiliaryMatrix")

        call zero(D_m)

        call allocate(diff_rhs, tfield%mesh, name=trim(field_name)//"DiffusionRHS")

        call zero(diff_rhs)

      end if

      ! allocate the rhs of the equation

      call allocate(rhs, tfield%mesh, name=trim(field_name)//"RHS")

      ! are we including a porosity coefficient on the time term?

      if (have_option(trim(complete_field_path(tfield%option_path))//'/porosity')) then

         include_porosity = .true.

         ! get the porosity theta averaged field - this will allocate it

         call form_porosity_theta(porosity_theta, state(1), &

             &option_path = trim(complete_field_path(tfield%option_path))//'/porosity')       

         call allocate(t_cvmass_with_porosity, tfield%mesh, name="CVMassWithPorosity")

         call compute_cv_mass(x, t_cvmass_with_porosity, porosity_theta)

         ewrite_minmax(t_cvmass_with_porosity)

         call deallocate(porosity_theta)

      else

         include_porosity = .false.

      end if

      ! find the cv mass that is used for the absorption and source terms

      t_abs_src_cvmass => get_cv_mass(state, tfield%mesh)

      ewrite_minmax(t_abs_src_cvmass)

      ! find the cv mass that is used for the time term derivative

      if (include_porosity) then

         t_cvmass => t_cvmass_with_porosity

      else

         t_cvmass => t_abs_src_cvmass      

      end if

      ewrite_minmax(t_cvmass)

      move_mesh = have_option("/mesh_adaptivity/mesh_movement")

      if(move_mesh) then

        if(.not.include_advection) then

          FLExit("Moving the mesh but not including advection is not possible yet.")

        end if

        if (include_porosity) then

           FLExit("Moving mesh not set up to work when including porosity")

        end if

        ewrite(2,*) "Moving mesh."

        x_old=>extract_vector_field(state(1), "OldCoordinate")

        x_new=>extract_vector_field(state(1), "IteratedCoordinate")

        call allocate(t_cvmass_old, tfield%mesh, name=trim(field_name)//"OldCVMass")

        call allocate(t_cvmass_new, tfield%mesh, name=trim(field_name)//"NewCVMass")

        call compute_cv_mass(x_old, t_cvmass_old)

        call compute_cv_mass(x_new, t_cvmass_new)

        ewrite_minmax(t_cvmass_old)

        ewrite_minmax(t_cvmass_new)

        ug=>extract_vector_field(state(1), "GridVelocity")

        ewrite_minmax(ug)

        ug_cvshape=make_cv_element_shape(cvfaces, ug%mesh%shape)

        ug_cvbdyshape=make_cvbdy_element_shape(cvfaces, ug%mesh%faces%shape)

      else

        ewrite(2,*) "Not moving mesh."

        ug_cvshape=u_cvshape

        ug_cvbdyshape=u_cvbdyshape

        call incref(ug_cvshape)

        call incref(ug_cvbdyshape)

      end if

      if(include_diffusion.and.(tfield_options%diffusionscheme==CV_DIFFUSION_BASSIREBAY)) then

        if(.not.(tfield%mesh==diffusivity%mesh)) then

          q_cvmass => get_cv_mass(state, diffusivity%mesh)

        else

          q_cvmass => t_cvmass

        end if

      else

        q_cvmass => t_cvmass

      end if

      ! allocate a field to store the locally iterated values in

      call allocate(advit_tfield, tfield%mesh, name="AdvIterated"//trim(field_name))

      ! allocate a field to use as the local old field for subcycling

      call allocate(l_old_tfield, tfield%mesh, name="LocalOld"//trim(field_name))

      ! allocate a field to store the change between the old and new values

      call allocate(delta_tfield, tfield%mesh, name="Delta_"//trim(field_name))

      call zero(delta_tfield) ! Impose zero initial guess.

      ! Ensure delta_tfield inherits options from tfield for solver

      delta_tfield%option_path = option_path

      ! find out how many iterations we'll be doing

      call get_option(trim(option_path)//"/prognostic/temporal_discretisation&

                      &/control_volumes/number_advection_iterations", &

                      adv_iterations, default=1)

      call get_option(trim(option_path)//"/prognostic/temporal_discretisation&

                      &/control_volumes/number_advection_iterations/tolerance", &

                      adv_tolerance, default=0.0)

      sub_dt=dt  ! just in case I don't initialise this somehow

      ! are we subcycling?

      no_subcycles = 1

      call get_option(trim(option_path)//"/prognostic/temporal_discretisation&

                      &/control_volumes/number_advection_subcycles", &

                      no_subcycles, stat=stat)

      if(stat/=0) then

        ! have not specified a number of subcycles but perhaps we're using a 

        ! courant number definition?

        call get_option(trim(option_path)//"/prognostic/temporal_discretisation&

                        &/control_volumes/maximum_courant_number_per_subcycle", &

                        max_sub_cfl, stat=stat)

        if(stat==0) then

          max_cfl = maxval(cfl_no%val)

          call allmax(max_cfl)

          ! yes, we're subcycling

          ! we should have already calculated the courant number (or aborted in the attempt)

          no_subcycles=ceiling(max_cfl/max_sub_cfl)

          if(no_subcycles>1) then

            sub_dt=dt/real(no_subcycles)

            call scale(cfl_no, 1.0/real(no_subcycles))

          end if

        else

          ! no, we're not subcycling

          no_subcycles=1

          sub_dt = dt

        end if

      else

        if(no_subcycles>1) then

          sub_dt=dt/real(no_subcycles)

          call scale(cfl_no, 1.0/real(no_subcycles))

        end if

      end if

      ! when subcycling we're going to need to be starting each subcycle from the

      ! "new" old value but I don't want to screw with old code by updating the actual

      ! global timestep old value so lets create a copy now and update it instead

      call set(l_old_tfield, oldtfield)

      ewrite(2,*) 'entering subcycling_loop', no_subcycles

      ! subcycling loop

      subcycling_loop: do sub = 1, no_subcycles

        ! advection iteration loop

        advection_iteration_loop: do adv_it = 1, adv_iterations

          ! construct the advection matrix if this is the first advection iteration

          ! and the first subcycle and we're not explicit and we're including advection

          assemble_advection_matrix=(adv_it==1).and.(sub==1).and.(.not.explicit).and.include_advection

          ! construct the diffusion matrix and rhs if we're including diffusion and

          ! if we're on the first advection iteration and the first subcycle

          assemble_diffusion=(adv_it==1).and.(sub==1).and.include_diffusion

          ! record the value of tfield since the previous iteration

          call set(advit_tfield, tfield)

          if(include_advection.or.assemble_diffusion) then

            ! if advection is being included or we need to assemble a

            ! diffusion matrix then assemble A_m, D_m and rhs here

            call assemble_advectiondiffusion_m_cv(A_m, rhs, D_m, diff_rhs, &

                                        tfield, l_old_tfield, tfield_options, &

                                        tdensity, oldtdensity, tdensity_options, &

                                        cvfaces, x_cvshape, x_cvbdyshape, &

                                        u_cvshape, u_cvbdyshape, t_cvshape, &

                                        ug_cvshape, ug_cvbdyshape, &

                                        x_cvshape_full, x_cvbdyshape_full, &

                                        t_cvshape_full, t_cvbdyshape_full, &