reaction_diffusion_equation_routines.f90

Go to the documentation of this file.

MODULE reaction_diffusion_equation_routines

  USE base_routines
  USE basis_routines
  USE boundary_conditions_routines
  USE comp_environment
  USE constants
  USE control_loop_routines
  USE distributed_matrix_vector
  USE domain_mappings
  USE equations_routines
  USE equations_mapping_routines
  USE equations_matrices_routines
  USE equations_set_constants
  USE field_routines
  USE input_output
  USE iso_varying_string
  USE kinds
  USE matrix_vector
#ifndef NOMPIMOD
  USE mpi
#endif
  USE problem_constants
  USE strings
  USE solver_routines
  USE timer
  USE types
  
  USE reaction_diffusion_io_routines

#include "macros.h"  

  IMPLICIT NONE

  PRIVATE

#ifdef NOMPIMOD
#include "mpif.h"
#endif

  !Module parameters

  !Module types

  !Module variables

  !Interfaces

  PUBLIC reactiondiffusion_equationssetsetup

  PUBLIC reactiondiffusion_equationssetsolutionmethodset
  
  PUBLIC reactiondiffusion_equationssetspecificationset

  PUBLIC reactiondiffusion_finiteelementcalculate
  
  PUBLIC reaction_diffusion_pre_solve
  
  PUBLIC reaction_diffusion_equation_problem_setup

  PUBLIC reactiondiffusion_problemspecificationset

  PUBLIC reaction_diffusion_post_solve

  PUBLIC reaction_diffusion_control_loop_post_loop


CONTAINS

  !
  !================================================================================================================================
  !

  SUBROUTINE reactiondiffusion_equationssetsetup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: component_idx,DIMENSION_MULTIPLIER,GEOMETRIC_COMPONENT_NUMBER,GEOMETRIC_SCALING_TYPE, &
      & NUMBER_OF_DIMENSIONS,NUMBER_OF_MATERIALS_COMPONENTS,GEOMETRIC_MESH_COMPONENT
    TYPE(decomposition_type), POINTER :: GEOMETRIC_DECOMPOSITION
    TYPE(equations_type), POINTER :: EQUATIONS
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_set_materials_type), POINTER :: EQUATIONS_MATERIALS
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("ReactionDiffusion_EquationsSetSetup",err,error,*999)

    NULLIFY(equations)
    NULLIFY(equations_mapping)
    NULLIFY(equations_matrices)
    NULLIFY(geometric_decomposition)

    IF(ASSOCIATED(equations_set)) THEN
      IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
        CALL flagerror("Equations set specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
        CALL flagerror("Equations set specification must have three entries for a reaction-diffusion type equations set.", &
          & err,error,*999)
      END IF
      SELECT CASE(equations_set_setup%SETUP_TYPE)
      CASE(equations_set_setup_initial_type)
        SELECT CASE(equations_set_setup%ACTION_TYPE)
        CASE(equations_set_setup_start_action)
          CALL reactiondiffusion_equationssetsolutionmethodset(equations_set, &
           & equations_set_fem_solution_method,err,error,*999)
        CASE(equations_set_setup_finish_action)
!!Todo: CHECK VALID SETUP
        CASE DEFAULT
          local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
            & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a reaction diffusion domain equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE(equations_set_setup_geometry_type)
        !\todo Check geometric dimension
      CASE(equations_set_setup_dependent_type)
        SELECT CASE(equations_set_setup%ACTION_TYPE)
        CASE(equations_set_setup_start_action)
          SELECT CASE(equations_set%SPECIFICATION(3))
          CASE(equations_set_cellml_reac_split_reac_diff_subtype, &
           & equations_set_constant_reac_diff_subtype)
            IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
              !Create the auto created dependent field
              CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_set%DEPENDENT% &
                & dependent_field,err,error,*999)
              CALL field_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,"Dependent Field",err,error,*999)
              CALL field_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_general_type,err,error,*999)
              CALL field_dependent_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_dependent_type,err,error,*999)
              CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
              CALL field_mesh_decomposition_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_decomposition, &
                & err,error,*999)
              CALL field_geometric_field_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,equations_set%GEOMETRY% &
                & geometric_field,err,error,*999)
              CALL field_number_of_variables_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,2,err,error,*999)
              CALL field_variable_types_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,[field_u_variable_type, &
                & field_deludeln_variable_type],err,error,*999)
              !CALL FIELD_VARIABLE_LABEL_SET
              CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
                & field_scalar_dimension_type,err,error,*999)
              CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
                & field_scalar_dimension_type,err,error,*999)
              CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
                & field_dp_type,err,error,*999)
              CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
                & field_dp_type,err,error,*999)
              CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
                & err,error,*999)
              CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
                & err,error,*999)
              !Default to the geometric interpolation setup
              CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type,1, &
                & geometric_mesh_component,err,error,*999)
              CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
                & geometric_mesh_component,err,error,*999)
              CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
                & geometric_mesh_component,err,error,*999)
              SELECT CASE(equations_set%SOLUTION_METHOD)
              CASE(equations_set_fem_solution_method)
                CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
                  & field_u_variable_type,1,field_node_based_interpolation,err,error,*999)
                CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
                  & field_deludeln_variable_type,1,field_node_based_interpolation,err,error,*999)
                !Default the scaling to the geometric field scaling
                CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
                CALL field_scaling_type_set(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_scaling_type,err,error,*999)
              CASE(equations_set_bem_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_fd_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_fv_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_gfem_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_gfv_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE DEFAULT
                local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
                  & " is invalid."
                CALL flagerror(local_error,err,error,*999)
              END SELECT
            ELSE
              !Check the user specified field
              CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
              CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
              CALL field_number_of_variables_check(equations_set_setup%FIELD,2,err,error,*999)
              CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type,field_deludeln_variable_type], &
                & err,error,*999)
              CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_scalar_dimension_type,err,error,*999)
              CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_scalar_dimension_type, &
                & err,error,*999)
              CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
              CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type,err,error,*999)
              CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,1,err,error,*999)
              CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type,1,err,error,*999)
              SELECT CASE(equations_set%SOLUTION_METHOD)
              CASE(equations_set_fem_solution_method)
                CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
                  & field_node_based_interpolation,err,error,*999)
                CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
                  & field_node_based_interpolation,err,error,*999)
              CASE(equations_set_bem_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_fd_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_fv_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_gfem_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_gfv_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE DEFAULT
                local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
                  & " is invalid."
                CALL flagerror(local_error,err,error,*999)
              END SELECT
            ENDIF
          CASE(equations_set_cellml_reac_no_split_reac_diff_subtype)
            CALL flagerror("Not implemented.",err,error,*999)
          CASE DEFAULT
            local_error="The equation set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
              & " is invalid for a reaction diffusion equation set class."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(equations_set_setup_finish_action)
          IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
            CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
          ENDIF
        CASE DEFAULT
          local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
            & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a reaction diffusion equation"
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE(equations_set_setup_materials_type)
        SELECT CASE(equations_set_setup%ACTION_TYPE)
        CASE(equations_set_setup_start_action)
          IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
            equations_materials=>equations_set%MATERIALS
            IF(ASSOCIATED(equations_materials)) THEN
              IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
                !Create the auto created materials field
                CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_materials% &
                  & materials_field,err,error,*999)
                CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
                CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type,err,error,*999)
                CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
                CALL field_mesh_decomposition_set_and_lock(equations_materials%MATERIALS_FIELD,geometric_decomposition, &
                  & err,error,*999)
                CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
                  & geometric_field,err,error,*999)
                CALL field_number_of_variables_set_and_lock(equations_materials%MATERIALS_FIELD,1,err,error,*999)
                CALL field_variable_types_set_and_lock(equations_materials%MATERIALS_FIELD,[field_u_variable_type], &
                  & err,error,*999)
                CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & field_vector_dimension_type,err,error,*999)
                CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & field_dp_type,err,error,*999)
                CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
                  & number_of_dimensions,err,error,*999)
                IF(equations_set%SPECIFICATION(3)==equations_set_cellml_reac_split_reac_diff_subtype) THEN
                  !Reaction Diffusion. Materials field components are 1 diffusion coeff for each dimension
                  !plus one for the storage coefficient in alpha(delC/delt) = Div(-kgradC)+cellmlRC
                  number_of_materials_components=number_of_dimensions+1
                  dimension_multiplier=1
                ELSEIF(equations_set%SPECIFICATION(3)==equations_set_constant_reac_diff_subtype) THEN
                  !Constant reaction + diffusion. Materials field has 1 diffuse coeff for each dimension
                  !plus one for the storage coefficient om  alpha(delC/delt) = Div(-kgradC)+const(x)_source
                  number_of_materials_components=number_of_dimensions+1
                  dimension_multiplier=1
                ENDIF
                !Set the number of materials components
                CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & number_of_materials_components,err,error,*999)
                !Default the first three materials components for diffusivity param to the first component geometric interpolation with const interpolation
                DO component_idx=1,number_of_dimensions
                  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
                    & 1,geometric_component_number,err,error,*999)
                  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                    & component_idx,geometric_component_number,err,error,*999)
                  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                    & component_idx,field_constant_interpolation,err,error,*999)
                ENDDO !components_idx
                !Default the storage co-efficient to the first geometric interpolation setup with constant interpolation
                IF(equations_set%SPECIFICATION(3)==equations_set_cellml_reac_split_reac_diff_subtype) THEN
                  component_idx=number_of_materials_components
                  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
                    & 1,geometric_component_number,err,error,*999)
                  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                    & component_idx,geometric_component_number,err,error,*999)
                  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                    & component_idx,field_constant_interpolation,err,error,*999)
                ENDIF
                !Default the field scaling to that of the geometric field
                CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
                CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
              ELSE
                !Check the user specified field
                CALL field_type_check(equations_set_setup%FIELD,field_material_type,err,error,*999)
                CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
                CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
                CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
                CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
                  & err,error,*999)
                CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
                CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
                  & number_of_dimensions,err,error,*999)
                IF(equations_set%SPECIFICATION(3)==equations_set_cellml_reac_split_reac_diff_subtype .OR. &
                 & equations_set%SPECIFICATION(3)==equations_set_constant_reac_diff_subtype) THEN
                  !Reaction Diffusion with cellml. Materials field components are 1 for storage coeff plus one for each dimension i.e., k
                  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,number_of_dimensions+1, &
                    & err,error,*999)
                ENDIF
              ENDIF
            ELSE
              CALL flagerror("Equations set materials is not associated.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
          ENDIF
        CASE(equations_set_setup_finish_action)
          equations_materials=>equations_set%MATERIALS
          IF(ASSOCIATED(equations_materials)) THEN
            IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
              !Finish creating the materials field
              CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
              !Set the default values for the materials field
              CALL field_number_of_components_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
                & number_of_dimensions,err,error,*999)
              IF(equations_set%SPECIFICATION(3)==equations_set_cellml_reac_split_reac_diff_subtype .OR. &
               & equations_set%SPECIFICATION(3)==equations_set_constant_reac_diff_subtype) THEN
                !Reaction Diffusion with cellml. Materials field components are 1 plus one for each dimension i.e.,storage coeff, and k.
                number_of_materials_components=number_of_dimensions+1
                dimension_multiplier=1
              ENDIF
              !set the diffusion coefficients to be 1.0
              DO component_idx=1,number_of_dimensions
                CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & field_values_set_type,component_idx,1.0_dp,err,error,*999)
              ENDDO !component_idx
               !Now set storage-coefficient
              component_idx=number_of_dimensions+1
              CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                & field_values_set_type,component_idx,1.0_dp,err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Equations set materials is not associated.",err,error,*999)
          ENDIF
        CASE DEFAULT
          local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
            & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a reaction diffusion equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE(equations_set_setup_source_type)
        SELECT CASE(equations_set_setup%ACTION_TYPE)
        CASE(equations_set_setup_start_action)
          IF(ASSOCIATED(equations_set%MATERIALS)) THEN
            IF(equations_set%MATERIALS%MATERIALS_FINISHED) THEN
              IF(ASSOCIATED(equations_set%SOURCE)) THEN
                IF(equations_set%SOURCE%SOURCE_FIELD_AUTO_CREATED) THEN
                  !Create the auto created source field
                  !Start field creation with name 'SOURCE_FIELD'
                  CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION, &
                    & equations_set%SOURCE%SOURCE_FIELD,err,error,*999)
                  !Create a general field
                  CALL field_type_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_general_type,err,error,*999)
                  !Label the field
                  CALL field_label_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,"Source Field",err,error,*999)
                  !Set the dependent type
                  CALL field_dependent_type_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_independent_type,err,error,*999)
                  !Set the field decomposition to be that of the geometric decomposition
                  CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
                  CALL field_mesh_decomposition_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,geometric_decomposition, &
                    & err,error,*999)
                  !Set the geometric field
                  CALL field_geometric_field_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,equations_set% & 
                    & geometry%GEOMETRIC_FIELD,err,error,*999)
                  !Set the field variables.
                  CALL field_number_of_variables_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,1,err,error,*999)
                  CALL field_variable_types_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,[field_u_variable_type],err,error,*999)
                  !Set the dimension
                  CALL field_dimension_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_u_variable_type, &
                    & field_scalar_dimension_type,err,error,*999)
                  !Set the data type
                  CALL field_data_type_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_u_variable_type, &
                    & field_dp_type,err,error,*999)
                  !Set the number of components to one
                  CALL field_number_of_components_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_u_variable_type,1, &
                    & err,error,*999)
                  !Get the geometric mesh component
                  CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, & 
                    & 1,geometric_mesh_component,err,error,*999)
                  !Default to the geometric interpolation setup
                  CALL field_component_mesh_component_set(equations_set%SOURCE%SOURCE_FIELD,field_u_variable_type,1, &
                    & geometric_mesh_component,err,error,*999)            
                  !Specify the interpolation to be same as geometric interpolation
                  SELECT CASE(equations_set%SOLUTION_METHOD)
                  CASE(equations_set_fem_solution_method)
                    CALL field_component_interpolation_set_and_lock(equations_set%SOURCE%SOURCE_FIELD,field_u_variable_type,1, &
                      & field_node_based_interpolation,err,error,*999)
                  CASE(equations_set_bem_solution_method)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(equations_set_fd_solution_method)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(equations_set_fv_solution_method)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(equations_set_gfem_solution_method)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(equations_set_gfv_solution_method)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE DEFAULT
                    local_error="The solution method of " &
                      & //trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// " is invalid."
                    CALL flagerror(local_error,err,error,*999)
                  END SELECT
                  !Set the scaling to be the same as the geometric field
                  CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type, &
                    & err,error,*999)
                  CALL field_scaling_type_set(equations_set%SOURCE%SOURCE_FIELD,geometric_scaling_type, &
                    & err,error,*999)
                ELSE
                  !Check the user specified field
                  CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
                  CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
                  CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
                  CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
                  CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_scalar_dimension_type, &
                    & err,error,*999)
                  CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
                  CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,1,err,error,*999)
                  
                  SELECT CASE(equations_set%SOLUTION_METHOD)
                  CASE(equations_set_fem_solution_method)
                    CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
                      & field_node_based_interpolation,err,error,*999)
                  CASE(equations_set_bem_solution_method)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(equations_set_fd_solution_method)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(equations_set_fv_solution_method)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(equations_set_gfem_solution_method)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(equations_set_gfv_solution_method)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE DEFAULT
                    local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD, &
                      &"*",err,error))//" is invalid."
                    CALL flagerror(local_error,err,error,*999)
                  END SELECT
                  
                ENDIF
              ELSE
                CALL flagerror("Equations set source is not associated.",err,error,*999)
              ENDIF
            ELSE
              CALL flagerror("Equations set materials field has not been finished.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Equations set materials is not associated.",err,error,*999)
          ENDIF
        CASE(equations_set_setup_finish_action)
          IF(equations_set%SOURCE%SOURCE_FIELD_AUTO_CREATED) THEN
            CALL field_create_finish(equations_set%SOURCE%SOURCE_FIELD,err,error,*999)
          ENDIF
        CASE DEFAULT
          local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
            & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a reaction diffusion equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE(equations_set_setup_analytic_type)
        SELECT CASE(equations_set_setup%ACTION_TYPE)
        CASE(equations_set_setup_start_action)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(equations_set_setup_finish_action)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE DEFAULT
          local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
            & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a reaction diffusion equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE(equations_set_setup_equations_type)
        SELECT CASE(equations_set_setup%ACTION_TYPE)
        CASE(equations_set_setup_start_action)
          IF(ASSOCIATED(equations_set%SOURCE)) THEN
            IF(equations_set%SOURCE%SOURCE_FINISHED) THEN
              !Create the equations
              CALL equations_create_start(equations_set,equations,err,error,*999)
              SELECT CASE(equations_set%SPECIFICATION(3))
              CASE(equations_set_cellml_reac_split_reac_diff_subtype)
                CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
              CASE(equations_set_constant_reac_diff_subtype)
                CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
              CASE(equations_set_cellml_reac_no_split_reac_diff_subtype)
                CALL equations_linearity_type_set(equations,equations_nonlinear,err,error,*999)
              CASE DEFAULT
                local_error="The equations matrices linearity set up of "// &
                  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
                CALL flagerror(local_error,err,error,*999)
              END SELECT
              CALL equations_time_dependence_type_set(equations,equations_first_order_dynamic,err,error,*999)
            ELSE
              CALL flagerror("Equations set source field has not been finished.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Equations set source is not associated.",err,error,*999)
          ENDIF
        CASE(equations_set_setup_finish_action)
          SELECT CASE(equations_set%SOLUTION_METHOD)
          CASE(equations_set_fem_solution_method)
            !Finish the creation of the equations
            CALL equations_set_equations_get(equations_set,equations,err,error,*999)
            CALL equations_create_finish(equations,err,error,*999)
            !Create the equations mapping.
            CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
            CALL equations_mapping_dynamic_matrices_set(equations_mapping,.true.,.true.,err,error,*999)
            CALL equations_mapping_dynamic_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
            CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type,err,error,*999)

            CALL equations_mapping_source_variable_type_set(equations_mapping,field_u_variable_type,err,error,*999)
            CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
            !Create the equations matrices
            CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
            !Set up matrix storage and structure
            IF(equations%LUMPING_TYPE==equations_lumped_matrices) THEN
              !Set up lumping
              CALL equations_matrices_dynamic_lumping_type_set(equations_matrices, &
                & [equations_matrix_unlumped,equations_matrix_lumped],err,error,*999)
              CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
                & [distributed_matrix_compressed_row_storage_type,distributed_matrix_diagonal_storage_type], &
                & err,error,*999)
              CALL equationsmatrices_dynamicstructuretypeset(equations_matrices, &
                [equations_matrix_fem_structure,equations_matrix_diagonal_structure],err,error,*999)
            ELSE
              SELECT CASE(equations%SPARSITY_TYPE)
              CASE(equations_matrices_full_matrices) 
                CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
                  & [distributed_matrix_block_storage_type,distributed_matrix_block_storage_type], &
                  & err,error,*999)
              CASE(equations_matrices_sparse_matrices)
                CALL equations_matrices_dynamic_storage_type_set(equations_matrices, &
                  & [distributed_matrix_compressed_row_storage_type,distributed_matrix_compressed_row_storage_type], &
                  & err,error,*999)
                CALL equationsmatrices_dynamicstructuretypeset(equations_matrices, &
                  [equations_matrix_fem_structure,equations_matrix_fem_structure],err,error,*999)                  
              CASE DEFAULT
                local_error="The equations matrices sparsity type of "// &
                  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
                CALL flagerror(local_error,err,error,*999)
              END SELECT
            ENDIF
            CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
          CASE(equations_set_bem_solution_method)
            CALL flagerror("Not implemented.",err,error,*999)
          CASE(equations_set_fd_solution_method)
            CALL flagerror("Not implemented.",err,error,*999)
          CASE(equations_set_fv_solution_method)
            CALL flagerror("Not implemented.",err,error,*999)
          CASE(equations_set_gfem_solution_method)
            CALL flagerror("Not implemented.",err,error,*999)
          CASE(equations_set_gfv_solution_method)
            CALL flagerror("Not implemented.",err,error,*999)
          CASE DEFAULT
            local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
              & " is invalid."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE DEFAULT
          local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
            & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a bioelectric domain equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE DEFAULT
        local_error="The setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
          & " is invalid for reaction diffusion equation."
        CALL flagerror(local_error,err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    ENDIF
    
    exits("ReactionDiffusion_EquationsSetSetup")
    RETURN
999 errorsexits("ReactionDiffusion_EquationsSetSetup",err,error)
    RETURN 1
    
  END SUBROUTINE reactiondiffusion_equationssetsetup

  !
  !================================================================================================================================
  !

  SUBROUTINE reactiondiffusion_equationssetsolutionmethodset(EQUATIONS_SET,SOLUTION_METHOD,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    INTEGER(INTG), INTENT(IN) :: SOLUTION_METHOD
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("ReactionDiffusion_EquationsSetSolutionMethodSet",err,error,*999)
    
    IF(ASSOCIATED(equations_set)) THEN
      IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
        CALL flagerror("Equations set specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
        CALL flagerror("Equations set specification must have three entries for a reaction-diffusion type equations set.", &
          & err,error,*999)
      END IF
      SELECT CASE(equations_set%SPECIFICATION(3))
      CASE(equations_set_cellml_reac_split_reac_diff_subtype, &
         & equations_set_cellml_reac_no_split_reac_diff_subtype, &
         & equations_set_constant_reac_diff_subtype)        
        SELECT CASE(solution_method)
        CASE(equations_set_fem_solution_method)
          equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
        CASE(equations_set_bem_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(equations_set_fd_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(equations_set_fv_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(equations_set_gfem_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(equations_set_gfv_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE DEFAULT
          local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE DEFAULT
        local_error="Equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
          & " is not valid for a reaction diffusion equation type of classical equations set class."
        CALL flagerror(local_error,err,error,*999)
      END SELECT

    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    ENDIF
       
    exits("ReactionDiffusion_EquationsSetSolutionMethodSet")
    RETURN
999 errors("ReactionDiffusion_EquationsSetSolutionMethodSet",err,error)
    exits("ReactionDiffusion_EquationsSetSolutionMethodSet")
    RETURN 1
    
  END SUBROUTINE reactiondiffusion_equationssetsolutionmethodset

  !
  !================================================================================================================================ 
  !

  SUBROUTINE reactiondiffusion_equationssetspecificationset(equationsSet,specification,err,error,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: equationsSet
    INTEGER(INTG), INTENT(IN) :: specification(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(varying_string) :: localError
    INTEGER(INTG) :: subtype

    enters("ReactionDiffusion_EquationsSetSpecificationSet",err,error,*999)

    IF(ASSOCIATED(equationsset)) THEN
      IF(SIZE(specification,1)>3) THEN
        CALL flagerror("Equations set specification must have 3 entries for a reaction-diffusion type equations set.", &
          & err,error,*999)
      END IF
      subtype=specification(3)
      SELECT CASE(subtype)
      CASE(equations_set_cellml_reac_split_reac_diff_subtype)
        !ok
      CASE(equations_set_cellml_reac_no_split_reac_diff_subtype)
        CALL flagerror("Not implemented.",err,error,*999)
      CASE(equations_set_constant_reac_diff_subtype)
        !ok
      CASE DEFAULT
        localerror="The specified equations set subtype of "//trim(numbertovstring(subtype,"*",err,error))// &
          & " is not valid for reaction diffusion equation type of a classical equations set class."
        CALL flagerror(localerror,err,error,*999)
      END SELECT
      !Set full specification
      IF(ALLOCATED(equationsset%specification)) THEN
        CALL flagerror("Equations set specification is already allocated.",err,error,*999)
      ELSE
        ALLOCATE(equationsset%specification(3),stat=err)
        IF(err/=0) CALL flagerror("Could not allocate equations set specification.",err,error,*999)
      END IF
      equationsset%specification(1:3)=[equations_set_classical_field_class,equations_set_reaction_diffusion_equation_type,subtype]
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    END IF

    exits("ReactionDiffusion_EquationsSetSpecificationSet")
    RETURN
999 errors("ReactionDiffusion_EquationsSetSpecificationSet",err,error)
    exits("ReactionDiffusion_EquationsSetSpecificationSet")
    RETURN 1
    
  END SUBROUTINE reactiondiffusion_equationssetspecificationset

  !
  !================================================================================================================================
  !
  SUBROUTINE reactiondiffusion_finiteelementcalculate(EQUATIONS_SET,ELEMENT_NUMBER,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) FIELD_VAR_TYPE,mh,mhs,ms,ng,nh,nhs,ni,nj,ns,component_idx
    LOGICAL :: USE_FIBRES
    REAL(DP) :: DIFFUSIVITY(3,3),DPHIDX(3,64),RWG,SUM,STORAGE_COEFFICIENT,C_PARAM
    TYPE(basis_type), POINTER :: DEPENDENT_BASIS,GEOMETRIC_BASIS,FIBRE_BASIS
    TYPE(equations_type), POINTER :: EQUATIONS
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_mapping_dynamic_type), POINTER :: DYNAMIC_MAPPING
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
    TYPE(equations_matrices_source_type), POINTER :: SOURCE_VECTOR
    TYPE(equations_matrix_type), POINTER :: DAMPING_MATRIX,STIFFNESS_MATRIX
    TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD,FIBRE_FIELD,MATERIALS_FIELD,SOURCE_FIELD
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE,GEOMETRIC_VARIABLE
    TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME
    
    enters("ReactionDiffusion_FiniteElementCalculate",err,error,*999)

    IF(ASSOCIATED(equations_set)) THEN
      IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
        CALL flagerror("Equations set specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
        CALL flagerror("Equations set specification must have three entries for a reaction-diffusion type equations set.", &
          & err,error,*999)
      END IF
      equations=>equations_set%EQUATIONS
      IF(ASSOCIATED(equations)) THEN
        
        dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
        geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
        materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
        IF(equations_set%SPECIFICATION(3)==equations_set_cellml_reac_no_split_reac_diff_subtype .OR. &
            & equations_set%SPECIFICATION(3)==equations_set_constant_reac_diff_subtype) THEN
           source_field=>equations%INTERPOLATION%SOURCE_FIELD
        ENDIF
        fibre_field=>equations%INTERPOLATION%FIBRE_FIELD
        use_fibres=ASSOCIATED(fibre_field)
        equations_mapping=>equations%EQUATIONS_MAPPING
        equations_matrices=>equations%EQUATIONS_MATRICES
        dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
          & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
        geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
          & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
        geometric_variable=>geometric_field%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR
        IF(use_fibres) fibre_basis=>fibre_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
          & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
        quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
        CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
          & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
        CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
          & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
        IF(use_fibres) CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations% &
          & interpolation%FIBRE_INTERP_PARAMETERS(field_u_variable_type)%PTR,err,error,*999)
        IF(equations_set%SPECIFICATION(3)==equations_set_cellml_reac_no_split_reac_diff_subtype .OR. &
            & equations_set%SPECIFICATION(3)==equations_set_constant_reac_diff_subtype) THEN
           CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
             & source_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
        ENDIF                  
        dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
        stiffness_matrix=>dynamic_matrices%MATRICES(1)%PTR
        damping_matrix=>dynamic_matrices%MATRICES(2)%PTR
        rhs_vector=>equations_matrices%RHS_VECTOR
        IF(equations_set%SPECIFICATION(3)==equations_set_cellml_reac_no_split_reac_diff_subtype .OR. &
            & equations_set%SPECIFICATION(3)==equations_set_constant_reac_diff_subtype) THEN
           source_vector=>equations_matrices%SOURCE_VECTOR
        ENDIF
        dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
        field_variable=>dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
        field_var_type=field_variable%VARIABLE_TYPE
        IF(stiffness_matrix%UPDATE_MATRIX.OR.damping_matrix%UPDATE_MATRIX.OR.rhs_vector%UPDATE_VECTOR) THEN
          !Loop over gauss points
          DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
            CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
              & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
            CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
              & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
            CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
              & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)
            IF(use_fibres) THEN
              CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
                & fibre_interp_point(field_u_variable_type)%PTR,err,error,*999)
              CALL field_interpolated_point_metrics_calculate(fibre_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
                & fibre_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
            ENDIF
            IF(equations_set%SPECIFICATION(3)==equations_set_cellml_reac_no_split_reac_diff_subtype .OR. &
            & equations_set%SPECIFICATION(3)==equations_set_constant_reac_diff_subtype) THEN
              CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
                & source_interp_point(field_u_variable_type)%PTR,err,error,*999)
            ENDIF
            !Calculate RWG.
            rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
              & quadrature_scheme%GAUSS_WEIGHTS(ng)
            !Calculate the diffusivity tensor
            diffusivity=0.0_dp
            IF(use_fibres) THEN
              !Calculate the diffusivity tensor in fibre coordinates
              CALL flagerror("Not implemented.",err,error,*999)
            ELSE
              !Use the diffusivity tensor in geometric coordinates
              DO nj=1,geometric_variable%NUMBER_OF_COMPONENTS !first three components of material field are the diffusivities
                diffusivity(nj,nj)=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(nj,1)
              ENDDO !nj
            ENDIF
            !Get the storage Coefficient, stored in the component after the diffusivities for each dimension
            component_idx=geometric_variable%NUMBER_OF_COMPONENTS+1
            storage_coefficient=equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(component_idx,1)
            !Compute basis dPhi/dx terms
            DO nj=1,geometric_variable%NUMBER_OF_COMPONENTS
              DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
                dphidx(nj,ms)=0.0_dp
                DO ni=1,dependent_basis%NUMBER_OF_XI
                  dphidx(nj,ms)=dphidx(nj,ms)+ &
                    & quadrature_scheme%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)* &
                    & equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%DXI_DX(ni,nj)
                ENDDO !ni
              ENDDO !ms
            ENDDO !nj            
            !Loop over field components
            mhs=0          
            DO mh=1,field_variable%NUMBER_OF_COMPONENTS
              !Loop over element rows
              DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
                mhs=mhs+1
                nhs=0
                !Loop over element columns
                DO nh=1,field_variable%NUMBER_OF_COMPONENTS
                  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
                    nhs=nhs+1
                    sum=0.0_dp
                    IF(stiffness_matrix%UPDATE_MATRIX) THEN
                      DO ni=1,geometric_variable%NUMBER_OF_COMPONENTS
                        DO nj=1,geometric_variable%NUMBER_OF_COMPONENTS
                          sum=sum+diffusivity(ni,nj)*dphidx(ni,mhs)*dphidx(nj,nhs)
                       ENDDO !nj
                      ENDDO !ni
                      stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+(sum*rwg)
                    ENDIF
                    IF(damping_matrix%UPDATE_MATRIX) THEN
                      damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+ &
                        & quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)* &
                        & quadrature_scheme%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)*storage_coefficient*rwg
                    ENDIF
                  ENDDO !ns
                ENDDO !nh

                IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
              ENDDO !ms
            ENDDO !mh
            IF(equations_set%SPECIFICATION(3)==equations_set_constant_reac_diff_subtype) THEN
              IF(source_vector%UPDATE_VECTOR) THEN
                c_param=equations%INTERPOLATION%SOURCE_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,no_part_deriv)
                mhs=0
                DO mh=1,field_variable%NUMBER_OF_COMPONENTS
                !Loop over element rows
                  DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
                    mhs=mhs+1
                    source_vector%ELEMENT_VECTOR%VECTOR(mhs)=source_vector%ELEMENT_VECTOR%VECTOR(mhs)+ &
                      & quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)*c_param*rwg
                  ENDDO
                ENDDO
              ENDIF
            ENDIF
            IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
          ENDDO !ng
        ENDIF
        !Scale factor adjustment
        IF(dependent_field%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
          CALL field_interpolationparametersscalefactorselementget(element_number,equations%INTERPOLATION% &
            & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
          mhs=0          
          DO mh=1,field_variable%NUMBER_OF_COMPONENTS
            !Loop over element rows
            DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
              mhs=mhs+1                    
              nhs=0
              IF(stiffness_matrix%UPDATE_MATRIX.OR.damping_matrix%UPDATE_MATRIX) THEN
                !Loop over element columns
                DO nh=1,field_variable%NUMBER_OF_COMPONENTS
                  DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
                    nhs=nhs+1
                    IF(stiffness_matrix%UPDATE_MATRIX) THEN
                      stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=stiffness_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
                        & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
                        & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
                    ENDIF
                    IF(damping_matrix%UPDATE_MATRIX) THEN
                      damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=damping_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
                        & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
                        & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
                    ENDIF
                  ENDDO !ns
                ENDDO !nh
              ENDIF
              IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)* &
                & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)
              IF(equations_set%SPECIFICATION(3)==equations_set_constant_reac_diff_subtype) THEN
                 IF(source_vector%UPDATE_VECTOR) source_vector%ELEMENT_VECTOR%VECTOR(mhs)= & 
                    & source_vector%ELEMENT_VECTOR%VECTOR(mhs)* &
                    & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)
              ENDIF
            ENDDO !ms
          ENDDO !mh
        ENDIF
      ELSE
        CALL flagerror("Equations set equations is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    ENDIF
       
    exits("ReactionDiffusion_FiniteElementCalculate")
    RETURN
999 errors("ReactionDiffusion_FiniteElementCalculate",err,error)
    exits("ReactionDiffusion_FiniteElementCalculate")
    RETURN 1
    
  END SUBROUTINE reactiondiffusion_finiteelementcalculate

  !
  !================================================================================================================================
  !
  SUBROUTINE reactiondiffusion_problemspecificationset(problem,problemSpecification,err,error,*)

    !Argument variables
    TYPE(problem_type), POINTER, INTENT(IN) :: problem
    INTEGER(INTG), INTENT(IN) :: problemSpecification(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(varying_string) :: localError
    INTEGER(INTG) :: problemSubtype

    enters("ReactionDiffusion_ProblemSpecificationSet",err,error,*999)

    IF(ASSOCIATED(problem)) THEN
      IF(SIZE(problemspecification,1)>=3) THEN
        problemsubtype=problemspecification(3)
        SELECT CASE(problemsubtype)
        CASE(problem_cellml_reac_integ_reac_diff_strang_split_subtype, &
            & problem_cellml_reac_eval_reac_diff_no_split_subtype, &
            & problem_constant_reac_diff_no_split_subtype)
          !ok
        CASE DEFAULT
          localerror="The specified problem subtype of "//trim(numbertovstring(problemsubtype,"*",err,error))// &
            & " is not valid for a reaction-diffusion problem type of a classical problem class."
          CALL flagerror(localerror,err,error,*999)
        END SELECT
        IF(ALLOCATED(problem%specification)) THEN
          CALL flagerror("Problem specification is already allocated.",err,error,*999)
        ELSE
          ALLOCATE(problem%specification(3),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate problem specification.",err,error,*999)
        END IF
        problem%specification(1:3)=[problem_classical_field_class,problem_reaction_diffusion_equation_type,problemsubtype]
      ELSE
        CALL flagerror("Reaction-diffusion problem specification must have >=3 entries.",err,error,*999)
      END IF
    ELSE
      CALL flagerror("Problem is not associated.",err,error,*999)
    END IF

    exits("ReactionDiffusion_ProblemSpecificationSet")
    RETURN
999 errors("ReactionDiffusion_ProblemSpecificationSet",err,error)
    exits("ReactionDiffusion_ProblemSpecificationSet")
    RETURN 1
    
  END SUBROUTINE reactiondiffusion_problemspecificationset

  !
  !================================================================================================================================
  !
  SUBROUTINE reaction_diffusion_equation_problem_setup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)

    !Argument variables
    TYPE(problem_type), POINTER :: PROBLEM
    TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(cellml_equations_type), POINTER :: CELLML_EQUATIONS
    TYPE(control_loop_type), POINTER :: CONTROL_LOOP,CONTROL_LOOP_ROOT
    TYPE(solver_type), POINTER :: SOLVER
    TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
    TYPE(solvers_type), POINTER :: SOLVERS
    TYPE(varying_string) :: LOCAL_ERROR

    NULLIFY(cellml_equations)
    NULLIFY(control_loop)
    NULLIFY(solver)
    NULLIFY(solvers)
    NULLIFY(solver_equations)
    
    enters("REACTION_DIFFUSION_EQUATION_PROBLEM_SETUP",err,error,*999)

    IF(ASSOCIATED(problem)) THEN
      IF(.NOT.ALLOCATED(problem%SPECIFICATION)) THEN
        CALL flagerror("Problem specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(problem%SPECIFICATION,1)<3) THEN
        CALL flagerror("Problem specification must have three entries for a reaction diffusion problem.",err,error,*999)
      END IF
      SELECT CASE(problem_setup%SETUP_TYPE)
      CASE(problem_setup_initial_type)
        SELECT CASE(problem_setup%ACTION_TYPE)
        CASE(problem_setup_start_action)
          !Do nothing????
        CASE(problem_setup_finish_action)
          !Do nothing????
        CASE DEFAULT
          local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
            & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a reaction diffusion equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE(problem_setup_control_type)
        SELECT CASE(problem_setup%ACTION_TYPE)
        CASE(problem_setup_start_action)
          !Set up a time control loop
          CALL control_loop_create_start(problem,control_loop,err,error,*999)
          CALL control_loop_type_set(control_loop,problem_control_time_loop_type,err,error,*999)
          CALL control_loop_label_set(control_loop,"Time Loop",err,error,*999)
        CASE(problem_setup_finish_action)
          !Finish the control loops
          control_loop_root=>problem%CONTROL_LOOP
          CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
          CALL control_loop_create_finish(control_loop,err,error,*999)            
        CASE DEFAULT
          local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
            & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a reaction-diffusion equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE(problem_setup_solvers_type)
        !Get the control loop
        control_loop_root=>problem%CONTROL_LOOP
        CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
        SELECT CASE(problem_setup%ACTION_TYPE)
        CASE(problem_setup_start_action)
          !Start the solvers creation
          CALL solvers_create_start(control_loop,solvers,err,error,*999)
          SELECT CASE(problem%SPECIFICATION(3))
          CASE(problem_cellml_reac_integ_reac_diff_strang_split_subtype)
            CALL solvers_number_set(solvers,3,err,error,*999)
            !Set the first solver to be a differential-algebraic equations solver
            NULLIFY(solver)
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            CALL solver_type_set(solver,solver_dae_type,err,error,*999)
            CALL solver_label_set(solver,"First ODE solver",err,error,*999)
            !Set solver defaults
            CALL solver_library_type_set(solver,solver_cmiss_library,err,error,*999)
            !Set the second solver to be a dynamic solver 
            NULLIFY(solver)
            CALL solvers_solver_get(solvers,2,solver,err,error,*999)
            CALL solver_type_set(solver,solver_dynamic_type,err,error,*999)
            CALL solver_dynamic_order_set(solver,solver_dynamic_first_order,err,error,*999)
            CALL solver_label_set(solver,"Parabolic solver",err,error,*999)
            !Set solver defaults
            CALL solver_dynamic_degree_set(solver,solver_dynamic_first_degree,err,error,*999)
            CALL solver_dynamic_scheme_set(solver,solver_dynamic_crank_nicolson_scheme,err,error,*999)
            CALL solver_library_type_set(solver,solver_cmiss_library,err,error,*999)
            !Set the third solver to be a differential-algebraic equations solver
            NULLIFY(solver)
            CALL solvers_solver_get(solvers,3,solver,err,error,*999)
            CALL solver_type_set(solver,solver_dae_type,err,error,*999)
            CALL solver_label_set(solver,"Second ODE solver",err,error,*999)
            !Set solver defaults
            CALL solver_library_type_set(solver,solver_cmiss_library,err,error,*999)
          CASE(problem_cellml_reac_eval_reac_diff_no_split_subtype)
            CALL solvers_number_set(solvers,2,err,error,*999)
            !Set the first solver to be a CELLML evaluator equations solver
            NULLIFY(solver)
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            CALL solver_type_set(solver,solver_cellml_evaluator_type,err,error,*999)
            CALL solver_label_set(solver,"Evaluator solver",err,error,*999)
            !Set solver defaults
            CALL solver_library_type_set(solver,solver_cmiss_library,err,error,*999)
            !Set the second solver to be a dynamic solver 
            NULLIFY(solver)
            CALL solvers_solver_get(solvers,2,solver,err,error,*999)
            CALL solver_type_set(solver,solver_dynamic_type,err,error,*999)
            CALL solver_dynamic_order_set(solver,solver_dynamic_first_order,err,error,*999)
            CALL solver_label_set(solver,"Parabolic solver",err,error,*999)
            !Set solver defaults
            CALL solver_dynamic_degree_set(solver,solver_dynamic_first_degree,err,error,*999)
            CALL solver_dynamic_scheme_set(solver,solver_dynamic_crank_nicolson_scheme,err,error,*999)
            CALL solver_library_type_set(solver,solver_cmiss_library,err,error,*999)
          CASE(problem_constant_reac_diff_no_split_subtype)
            CALL solvers_number_set(solvers,1,err,error,*999)
            !Set the solver to be a dynamic solver 
            NULLIFY(solver)
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            CALL solver_type_set(solver,solver_dynamic_type,err,error,*999)
            CALL solver_dynamic_order_set(solver,solver_dynamic_first_order,err,error,*999)
            CALL solver_label_set(solver,"Parabolic solver",err,error,*999)
            !Set solver defaults
            CALL solver_dynamic_degree_set(solver,solver_dynamic_first_degree,err,error,*999)
            CALL solver_dynamic_scheme_set(solver,solver_dynamic_crank_nicolson_scheme,err,error,*999)
            CALL solver_library_type_set(solver,solver_cmiss_library,err,error,*999)
          CASE DEFAULT
            local_error="The problem subtype of "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
              & " is invalid for a reaction-diffusion problem type of a classical problem class."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(problem_setup_finish_action)
          !Get the solvers
          CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
          !Finish the solvers creation
          CALL solvers_create_finish(solvers,err,error,*999)
        CASE DEFAULT
          local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
            & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a classical equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE(problem_setup_solver_equations_type)
        SELECT CASE(problem_setup%ACTION_TYPE)
        CASE(problem_setup_start_action)
          !Get the control loop
          control_loop_root=>problem%CONTROL_LOOP
          CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
          SELECT CASE(problem%SPECIFICATION(3))
          CASE(problem_cellml_reac_integ_reac_diff_strang_split_subtype)
            !Get the solver
            CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
            !Create the solver equations for the second (parabolic) solver
            CALL solvers_solver_get(solvers,2,solver,err,error,*999)
            CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
            CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
            CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_first_order_dynamic,err,error,*999)
            CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
          CASE(problem_cellml_reac_eval_reac_diff_no_split_subtype)
            !Get the solver
            CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
            !Create the solver equations for the parabolic solver
            CALL solvers_solver_get(solvers,2,solver,err,error,*999)
            CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
            CALL solver_equations_linearity_type_set(solver_equations,solver_equations_nonlinear,err,error,*999)
            CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_first_order_dynamic,err,error,*999)
            CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
          CASE(problem_constant_reac_diff_no_split_subtype)
            !Get the solver
            CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
            !Create the solver equations for the parabolic solver
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
            CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
            CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_first_order_dynamic,err,error,*999)
            CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup subtype of "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
              & " is invalid for a reaction-diffusion equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT

        CASE(problem_setup_finish_action)
          !Get the control loop
          control_loop_root=>problem%CONTROL_LOOP
          CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
          CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
          SELECT CASE(problem%SPECIFICATION(3))
          CASE(problem_cellml_reac_integ_reac_diff_strang_split_subtype)
            !Get the solver equations for the second (parabolic) solver
            CALL solvers_solver_get(solvers,2,solver,err,error,*999)
            CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
            !Finish the solver equations creation
            CALL solver_equations_create_finish(solver_equations,err,error,*999)             
          CASE(problem_cellml_reac_eval_reac_diff_no_split_subtype)
            !Get the solver equations for the parabolic solver
            CALL solvers_solver_get(solvers,2,solver,err,error,*999)
            CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
            !Finish the solver equations creation
            CALL solver_equations_create_finish(solver_equations,err,error,*999)             
          CASE(problem_constant_reac_diff_no_split_subtype)
            !Get the solver equations for thE solver
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
            !Finish the solver equations creation
            CALL solver_equations_create_finish(solver_equations,err,error,*999)             
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup subtype of "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
              & " is invalid for a reaction-diffusion equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE DEFAULT
          local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
            & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a reaction-diffusion equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE(problem_setup_cellml_equations_type)
        SELECT CASE(problem_setup%ACTION_TYPE)
        CASE(problem_setup_start_action)
          !Get the control loop
          control_loop_root=>problem%CONTROL_LOOP
          CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
          !Get the solver
          CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)

          IF(problem%SPECIFICATION(3)==problem_cellml_reac_integ_reac_diff_strang_split_subtype) THEN
            NULLIFY(solver)
            NULLIFY(cellml_equations)
            !Create the CellML equations for the first DAE solver
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            CALL cellml_equations_create_start(solver,cellml_equations,err,error,*999)
            !Create the CellML equations for the second DAE solver
            NULLIFY(solver)
            NULLIFY(cellml_equations)
            CALL solvers_solver_get(solvers,3,solver,err,error,*999)
            CALL cellml_equations_create_start(solver,cellml_equations,err,error,*999)
          ELSEIF(problem%SPECIFICATION(3)== &
           & problem_cellml_reac_eval_reac_diff_no_split_subtype) THEN
            !CREATE the CellML equations for the first evaluator solver
            NULLIFY(solver)
            NULLIFY(cellml_equations)
            !Create the CellML equations for the first cellml evaluator solver
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            CALL cellml_equations_create_start(solver,cellml_equations,err,error,*999)
          ENDIF
        CASE(problem_setup_finish_action)
          !Get the control loop
          control_loop_root=>problem%CONTROL_LOOP
          CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
          CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
          SELECT CASE(problem%SPECIFICATION(3))
          CASE(problem_cellml_reac_integ_reac_diff_strang_split_subtype)
            !Get the CellML equations for the first DAE solver
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            CALL solver_cellml_equations_get(solver,cellml_equations,err,error,*999)
            !Finish the CellML equations creation
            CALL cellml_equations_create_finish(cellml_equations,err,error,*999)
            !Get the CellML equations for the second DAE solver
            NULLIFY(solver)
            NULLIFY(cellml_equations)
            CALL solvers_solver_get(solvers,3,solver,err,error,*999)
            CALL solver_cellml_equations_get(solver,cellml_equations,err,error,*999)
            !Finish the CellML equations creation
            CALL cellml_equations_create_finish(cellml_equations,err,error,*999)
          CASE(problem_cellml_reac_eval_reac_diff_no_split_subtype)
            !Get the CellML equations for the first evaluator solver
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            CALL solver_cellml_equations_get(solver,cellml_equations,err,error,*999)
            !Finish the CellML equations creation
            CALL cellml_equations_create_finish(cellml_equations,err,error,*999)
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
              & " is invalid for reaction-diffusion equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE DEFAULT
          local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
            & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for reaction-diffusion equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE DEFAULT
        local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
          & " is invalid for areaction-diffusion equation."
        CALL flagerror(local_error,err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Problem is not associated.",err,error,*999)
    ENDIF
      
    exits("REACTION_DIFFUSION_EQUATION_PROBLEM_SETUP")
    RETURN
999 errorsexits("REACTION_DIFFUSION_EQUATION_PROBLEM_SETUP",err,error)
    RETURN 1
  END SUBROUTINE reaction_diffusion_equation_problem_setup
  
  !
  !================================================================================================================================
  !
  SUBROUTINE reaction_diffusion_pre_solve(SOLVER,ERR,ERROR,*)

    !Argument variables
    TYPE(solver_type), POINTER :: SOLVER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    REAL(DP) :: CURRENT_TIME,TIME_INCREMENT
    TYPE(control_loop_type), POINTER :: CONTROL_LOOP
    TYPE(problem_type), POINTER :: PROBLEM
    TYPE(solvers_type), POINTER :: SOLVERS
    TYPE(varying_string) :: LOCAL_ERROR

    enters("REACTION_DIFFUSION_PRE_SOLVE",err,error,*999)

    IF(ASSOCIATED(solver)) THEN
      solvers=>solver%SOLVERS
      IF(ASSOCIATED(solvers)) THEN
        control_loop=>solvers%CONTROL_LOOP
        IF(ASSOCIATED(control_loop)) THEN
          CALL control_loop_current_times_get(control_loop,current_time,time_increment,err,error,*999)
          problem=>control_loop%PROBLEM
          IF(ASSOCIATED(problem)) THEN
            SELECT CASE(problem%SPECIFICATION(3))
            CASE(problem_cellml_reac_integ_reac_diff_strang_split_subtype)
              SELECT CASE(solver%GLOBAL_NUMBER)
              CASE(1)
                CALL solver_dae_times_set(solver,current_time,current_time+time_increment/2.0_dp,err,error,*999)
              CASE(2)
                !Do nothing
              CASE(3)
                CALL solver_dae_times_set(solver,current_time,current_time+time_increment/2.0_dp, &
                  & err,error,*999)
              CASE DEFAULT
                local_error="The solver global number of "//trim(number_to_vstring(solver%GLOBAL_NUMBER,"*",err,error))// &
                  & " is invalid for a Strang split reaction-diffusion problem."
                CALL flagerror(local_error,err,error,*999)
              END SELECT
            CASE(problem_constant_reac_diff_no_split_subtype)
              !No splitting, therefore entire problem is solved as a dynamic one, with 1 solver nothing to do.
            CASE(problem_cellml_reac_eval_reac_diff_no_split_subtype)
              !No splitting, therefore entire problem is solved as a dynamic one, with 1 solver nothing to do.
            CASE DEFAULT
              local_error="The problem subtype of "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
                & " is invalid for a reaction-diffusion problem type."
              CALL flagerror(local_error,err,error,*999)
            END SELECT
          ELSE
            CALL flagerror("Control loop problem is not associated.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Solvers control loop is not associated.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Solver solvers is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Solver is not associated.",err,error,*999)
    ENDIF
       
    exits("REACTION_DIFFUSION_PRE_SOLVE")
    RETURN
999 errorsexits("REACTION_DIFFUSION_PRE_SOLVE",err,error)
    RETURN 1
    
  END SUBROUTINE reaction_diffusion_pre_solve

  !
  !================================================================================================================================
  !

  SUBROUTINE reaction_diffusion_post_solve(CONTROL_LOOP,SOLVER,ERR,ERROR,*)

    !Argument variables
    TYPE(control_loop_type), POINTER :: CONTROL_LOOP
    TYPE(solver_type), POINTER :: SOLVER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(solvers_type), POINTER :: SOLVERS
    TYPE(solver_type), POINTER :: PDE_SOLVER
    TYPE(varying_string) :: LOCAL_ERROR

    enters("REACTION_DIFFUSION_POST_SOLVE",err,error,*999)

    IF(ASSOCIATED(control_loop)) THEN
      IF(ASSOCIATED(solver)) THEN
        IF(ASSOCIATED(control_loop%PROBLEM)) THEN 
          IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
            CALL flagerror("Problem specification is not allocated.",err,error,*999)
          ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
            CALL flagerror("Problem specification must have three entries for a reaction diffusion problem.",err,error,*999)
          END IF
          SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
            CASE(problem_cellml_reac_integ_reac_diff_strang_split_subtype)
              SELECT CASE(solver%GLOBAL_NUMBER)
                CASE(1)

                !do nothing
                CASE(2)

                !do nothing
                CASE(3)
                  !OUTPUT SOLUTIONS AT EACH TIME STEP - should probably change this bit below to output 
                  !mesh solutions directly from the 3rd solver itself rather than by getting the 2nd solver.
                  !I just don't know how to work with cellml_equations to do this.
                  solvers=>solver%SOLVERS
                  NULLIFY(pde_solver)
                  CALL solvers_solver_get(solvers,2,pde_solver,err,error,*999)
                  CALL reaction_diffusion_post_solve_output_data(control_loop,pde_solver,err,error,*999)
                CASE DEFAULT
                  local_error="The solver global number of "//trim(number_to_vstring(solver%GLOBAL_NUMBER,"*",err,error))// &
                    & " is invalid for a Strang split reaction-diffusion problem."
                  CALL flagerror(local_error,err,error,*999)
              END SELECT
            CASE (problem_cellml_reac_eval_reac_diff_no_split_subtype)
              !do nothing - time output not implemented
            CASE (problem_constant_reac_diff_no_split_subtype)
              !OUTPUT SOLUTIONS AT TIME STEP
              CALL reaction_diffusion_post_solve_output_data(control_loop,solver,err,error,*999)
            CASE DEFAULT
              local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
                & " is not valid for a reaction diffusion type of a classical field problem class."
              CALL flagerror(local_error,err,error,*999)
          END SELECT
        ELSE
          CALL flagerror("Problem is not associated.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Solver is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Control loop is not associated.",err,error,*999)
    ENDIF

    exits("REACTION_DIFFUSION_POST_SOLVE")
    RETURN
999 errorsexits("REACTION_DIFFUSION_POST_SOLVE",err,error)
    RETURN 1
  END SUBROUTINE reaction_diffusion_post_solve
  !   
  !================================================================================================================================
  !

  SUBROUTINE reaction_diffusion_post_solve_output_data(CONTROL_LOOP,SOLVER,ERR,ERROR,*)

    !Argument variables
    TYPE(control_loop_type), POINTER :: CONTROL_LOOP
    TYPE(solver_type), POINTER :: SOLVER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
    TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    TYPE(varying_string) :: LOCAL_ERROR

    REAL(DP) :: CURRENT_TIME,TIME_INCREMENT
    INTEGER(INTG) :: EQUATIONS_SET_IDX,CURRENT_LOOP_ITERATION,OUTPUT_FREQUENCY,MAX_DIGITS
    INTEGER(INTG) :: myComputationalNodeNumber

    CHARACTER(30) :: FILE
    CHARACTER(30) :: OUTPUT_FILE
    
    CHARACTER(100) :: FMT, TEMP_FMT

    enters("REACTION_DIFFUSION_POST_SOLVE_OUTPUT_DATA",err,error,*999)

    IF(ASSOCIATED(control_loop)) THEN
      IF(ASSOCIATED(solver)) THEN
        IF(ASSOCIATED(control_loop%PROBLEM)) THEN
          IF(.NOT.ALLOCATED(control_loop%PROBLEM%SPECIFICATION)) THEN
            CALL flagerror("Problem specification is not allocated.",err,error,*999)
          ELSE IF(SIZE(control_loop%PROBLEM%SPECIFICATION,1)<3) THEN
            CALL flagerror("Problem specification must have three entries for a reaction diffusion problem.",err,error,*999)
          END IF
          SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
            CASE(problem_cellml_reac_integ_reac_diff_strang_split_subtype, &
              & problem_constant_reac_diff_no_split_subtype)
              CALL control_loop_current_times_get(control_loop,current_time,time_increment,err,error,*999)
              solver_equations=>solver%SOLVER_EQUATIONS
              IF(ASSOCIATED(solver_equations)) THEN
                solver_mapping=>solver_equations%SOLVER_MAPPING
                IF(ASSOCIATED(solver_mapping)) THEN
                  !Make sure the equations sets are up to date
                  DO equations_set_idx=1,solver_mapping%NUMBER_OF_EQUATIONS_SETS
                    equations_set=>solver_mapping%EQUATIONS_SETS(equations_set_idx)%PTR
                    current_loop_iteration=control_loop%TIME_LOOP%ITERATION_NUMBER
                    output_frequency=control_loop%TIME_LOOP%OUTPUT_NUMBER
                    mycomputationalnodenumber = computational_node_number_get(err,error)
                    max_digits=floor(log10((control_loop%TIME_LOOP%STOP_TIME-control_loop%TIME_LOOP%START_TIME)/ &
                      & control_loop%TIME_LOOP%TIME_INCREMENT))+1
                    IF(output_frequency>0) THEN
                      IF(mod(current_loop_iteration,output_frequency)==0) THEN
                        IF(control_loop%TIME_LOOP%CURRENT_TIME<=control_loop%TIME_LOOP%STOP_TIME) THEN
                          IF(solver_mapping%NUMBER_OF_EQUATIONS_SETS.EQ.1) THEN
                            WRITE(temp_fmt,'("I",I0,".",I0)') max_digits,max_digits
                            !100 FORMAT 
                            fmt = trim(temp_fmt)
                            WRITE(temp_fmt,'(A2,A38,A20,A2)') "(", '"TIME_STEP_SPEC_1.part",I2.2,".",',fmt,")"
                            fmt = trim(temp_fmt)
                            WRITE(output_file,fmt) &
                              & mycomputationalnodenumber,current_loop_iteration
                          ELSE
                            WRITE(temp_fmt,'("I",I0,".",I0)') max_digits,max_digits
                            !200 FORMAT 
                            fmt = trim(temp_fmt)
                            WRITE(temp_fmt,'(A2,A38,A20,A2)') "(", '"TIME_STEP_SPEC_",I0,".part",I2.2,".",',fmt,")"
                            fmt = trim(temp_fmt)
                            WRITE(output_file,fmt) &
                              & equations_set_idx, mycomputationalnodenumber,current_loop_iteration
                          ENDIF
                          WRITE(*,*) output_file
                          file=trim(output_file)
                          CALL write_string(general_output_type,"...",err,error,*999)
                          CALL write_string(general_output_type,"Now export fields... ",err,error,*999)
                          CALL reaction_diffusion_io_write_cmgui(equations_set%REGION,equations_set%GLOBAL_NUMBER,file, &
                            & err,error,*999)
                        ENDIF
                      ENDIF 
                    ENDIF
                  ENDDO
                ENDIF
              ENDIF
            CASE(problem_cellml_reac_eval_reac_diff_no_split_subtype)
              ! do nothing ???
              CALL flagerror("Not implemented.",err,error,*999)
            CASE DEFAULT
              local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
                & " is not valid for an advection-diffusion equation type of a classical field problem class."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        ELSE
          CALL flagerror("Problem is not associated.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Solver is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Control loop is not associated.",err,error,*999)
    ENDIF

    exits("REACTION_DIFFUSION_POST_SOLVE_OUTPUT_DATA")
    RETURN
999 errorsexits("REACTION_DIFFUSION_POST_SOLVE_OUTPUT_DATA",err,error)
    RETURN 1
    
  END SUBROUTINE reaction_diffusion_post_solve_output_data
  !
  !================================================================================================================================
  !

  SUBROUTINE reaction_diffusion_control_loop_post_loop(CONTROL_LOOP,ERR,ERROR,*)

    !Argument variables
    TYPE(control_loop_type), POINTER :: CONTROL_LOOP
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(problem_type), POINTER :: PROBLEM
    TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
    TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
    TYPE(equations_type), POINTER :: EQUATIONS
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    TYPE(solver_type), POINTER :: SOLVER
    TYPE(solvers_type), POINTER :: SOLVERS
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    TYPE(equations_matrix_type), POINTER :: DAMPING_MATRIX,STIFFNESS_MATRIX
    
    enters("REACTION_DIFFUSION_CONTROL_LOOP_POST_LOOP",err,error,*999)
    
    NULLIFY(solver)
    IF(ASSOCIATED(control_loop)) THEN
      problem=>control_loop%PROBLEM
      IF(ASSOCIATED(problem)) THEN
        SELECT CASE(problem%SPECIFICATION(3))
        CASE(problem_cellml_reac_integ_reac_diff_strang_split_subtype)
          solvers=>control_loop%SOLVERS
          IF(ASSOCIATED(solvers)) THEN
            CALL solvers_solver_get(solvers,2,solver,err,error,*999)
            solver_equations=>solver%SOLVER_EQUATIONS
            IF(ASSOCIATED(solver_equations)) THEN
              solver_mapping=>solver_equations%SOLVER_MAPPING
              IF(ASSOCIATED(solver_mapping)) THEN
                equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
                IF(ASSOCIATED(equations_set)) THEN
                  equations=>equations_set%EQUATIONS
                  IF(ASSOCIATED(equations)) THEN
                    equations_matrices=>equations%EQUATIONS_MATRICES
                    dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
                    stiffness_matrix=>dynamic_matrices%MATRICES(1)%PTR
                    damping_matrix=>dynamic_matrices%MATRICES(2)%PTR
                    stiffness_matrix%UPDATE_MATRIX = .false.
                    damping_matrix%UPDATE_MATRIX = .false.
                  ELSE
                    CALL flagerror("Equations not associated.",err,error,*999)
                  ENDIF
                ELSE
                  CALL flagerror("Equations Set not associated.",err,error,*999)
                ENDIF
      
              ELSE
                CALL flagerror("Solver Mapping not associated.",err,error,*999)
              ENDIF
            ELSE
              CALL flagerror("Solver Equations not associated.", err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Solvers is not associated.", err,error,*999)
          ENDIF


        CASE DEFAULT
          !do nothing
        END SELECT
      ELSE
        CALL flagerror("Problem is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Control Loop is not associated.",err,error,*999)
    ENDIF
    exits("REACTION_DIFFUSION_CONTROL_LOOP_POST_LOOP")
    RETURN
999 errorsexits("REACTION_DIFFUSION_CONTROL_LOOP_POST_LOOP",err,error)
    RETURN 1
  END SUBROUTINE reaction_diffusion_control_loop_post_loop
END MODULE reaction_diffusion_equation_routines