equations_matrices_routines.f90

Go to the documentation of this file.

MODULE equations_matrices_routines

  USE base_routines
  USE distributed_matrix_vector
  USE equations_set_constants
  USE field_routines
  USE iso_varying_string
  USE kinds
  USE lists
  USE matrix_vector
  USE strings
  USE types
  USE linkedlist_routines

#include "macros.h"  

  IMPLICIT NONE

  PRIVATE

  !Module parameters

  INTEGER(INTG), PARAMETER :: equations_matrix_no_structure=1
  INTEGER(INTG), PARAMETER :: equations_matrix_fem_structure=2
  INTEGER(INTG), PARAMETER :: equations_matrix_diagonal_structure=3
  INTEGER(INTG), PARAMETER :: equations_matrix_nodal_structure=4


  INTEGER(INTG), PARAMETER :: equations_matrix_unlumped=1
  INTEGER(INTG), PARAMETER :: equations_matrix_lumped=2
 
  INTEGER(INTG), PARAMETER :: equations_matrices_sparse_matrices=1
  INTEGER(INTG), PARAMETER :: equations_matrices_full_matrices=2

  INTEGER(INTG), PARAMETER :: equations_matrices_all=1
  INTEGER(INTG), PARAMETER :: equations_matrices_dynamic_only=2
  INTEGER(INTG), PARAMETER :: equations_matrices_linear_only=3
  INTEGER(INTG), PARAMETER :: equations_matrices_nonlinear_only=4
  INTEGER(INTG), PARAMETER :: equations_matrices_jacobian_only=5
  INTEGER(INTG), PARAMETER :: equations_matrices_residual_only=6
  INTEGER(INTG), PARAMETER :: equations_matrices_rhs_only=7
  INTEGER(INTG), PARAMETER :: equations_matrices_source_only=8
  INTEGER(INTG), PARAMETER :: equations_matrices_rhs_residual_only=9
  INTEGER(INTG), PARAMETER :: equations_matrices_rhs_source_only=10
  INTEGER(INTG), PARAMETER :: equations_matrices_residual_source_only=11
  INTEGER(INTG), PARAMETER :: equations_matrices_vectors_only=12

  INTEGER(INTG), PARAMETER :: equations_jacobian_finite_difference_calculated=1
  INTEGER(INTG), PARAMETER :: equations_jacobian_analytic_calculated=2

  !Module types

  !Module variables

  !Interfaces

  INTERFACE equationsmatrices_createfinish
    MODULE PROCEDURE equations_matrices_create_finish
  END INTERFACE equationsmatrices_createfinish

  INTERFACE equationsmatrices_createstart
    MODULE PROCEDURE equations_matrices_create_start
  END INTERFACE equationsmatrices_createstart

  INTERFACE equationsmatrices_destroy
    MODULE PROCEDURE equations_matrices_destroy
  END INTERFACE equationsmatrices_destroy

  INTERFACE equationsmatrices_dynamiclumpingtypeset
    MODULE PROCEDURE equations_matrices_dynamic_lumping_type_set
  END INTERFACE equationsmatrices_dynamiclumpingtypeset

  INTERFACE equationsmatrices_dynamicstoragetypeset
    MODULE PROCEDURE equations_matrices_dynamic_storage_type_set
  END INTERFACE equationsmatrices_dynamicstoragetypeset

  INTERFACE equationsmatrices_elementinitialise
    MODULE PROCEDURE equations_matrices_element_initialise
  END INTERFACE equationsmatrices_elementinitialise

  INTERFACE equationsmatrices_elementfinalise
    MODULE PROCEDURE equations_matrices_element_finalise
  END INTERFACE equationsmatrices_elementfinalise

  INTERFACE equationsmatrices_elementadd
    MODULE PROCEDURE equations_matrices_element_add
  END INTERFACE equationsmatrices_elementadd
  
  INTERFACE equationsmatrices_jacobianelementadd
    MODULE PROCEDURE equations_matrices_jacobian_element_add
  END INTERFACE equationsmatrices_jacobianelementadd
  
  INTERFACE equationsmatrices_elementcalculate
    MODULE PROCEDURE equations_matrices_element_calculate
  END INTERFACE equationsmatrices_elementcalculate

  INTERFACE equationsmatrices_valuesinitialise
    MODULE PROCEDURE equations_matrices_values_initialise
  END INTERFACE equationsmatrices_valuesinitialise

  INTERFACE equationsmatrices_elementmatrixfinalise
    MODULE PROCEDURE equations_matrices_element_matrix_finalise
  END INTERFACE equationsmatrices_elementmatrixfinalise

  INTERFACE equationsmatrices_elementmatrixcalculate
    MODULE PROCEDURE equations_matrices_element_matrix_calculate
  END INTERFACE equationsmatrices_elementmatrixcalculate

  INTERFACE equationsmatrices_elementmatrixsetup
    MODULE PROCEDURE equations_matrices_element_matrix_setup
  END INTERFACE equationsmatrices_elementmatrixsetup

  INTERFACE equationsmatrices_elementvectorfinalise
    MODULE PROCEDURE equations_matrices_element_vector_finalise
  END INTERFACE equationsmatrices_elementvectorfinalise

  INTERFACE equationsmatrices_elementvectorcalculate
    MODULE PROCEDURE equations_matrices_element_vector_calculate
  END INTERFACE equationsmatrices_elementvectorcalculate

  INTERFACE equationsmatrices_elementvectorsetup
    MODULE PROCEDURE equations_matrices_element_vector_setup
  END INTERFACE equationsmatrices_elementvectorsetup

  INTERFACE equationsmatrices_linearstoragetypeset
    MODULE PROCEDURE equations_matrices_linear_storage_type_set
  END INTERFACE equationsmatrices_linearstoragetypeset

  INTERFACE equationsmatrices_nonlinearstoragetypeset
    MODULE PROCEDURE equationsmatrices_nonlinearstoragetypeset0
    MODULE PROCEDURE equationsmatrices_nonlinearstoragetypeset1
  END INTERFACE equationsmatrices_nonlinearstoragetypeset

  INTERFACE equationsmatrices_nonlinearstructuretypeset
    MODULE PROCEDURE equationsmatrices_nonlinearstructuretypeset0
    MODULE PROCEDURE equationsmatrices_nonlinearstructuretypeset1
  END INTERFACE equationsmatrices_nonlinearstructuretypeset

  INTERFACE equationsmatrices_output
    MODULE PROCEDURE equations_matrices_output
  END INTERFACE equationsmatrices_output

  INTERFACE equationsmatrices_jacobianoutput
    MODULE PROCEDURE equations_matrices_jacobian_output
  END INTERFACE equationsmatrices_jacobianoutput

  PUBLIC equations_matrix_no_structure,equations_matrix_fem_structure,equations_matrix_diagonal_structure, &
    & equations_matrix_nodal_structure

  PUBLIC equations_matrix_unlumped,equations_matrix_lumped

  PUBLIC equations_matrices_sparse_matrices,equations_matrices_full_matrices

  PUBLIC equations_matrices_all,equations_matrices_linear_only,equations_matrices_nonlinear_only,equations_matrices_jacobian_only, &
    & equations_matrices_residual_only,equations_matrices_rhs_only,equations_matrices_source_only, &
    & equations_matrices_rhs_residual_only,equations_matrices_rhs_source_only,equations_matrices_residual_source_only, &
    & equations_matrices_vectors_only
  
  PUBLIC equations_jacobian_finite_difference_calculated,equations_jacobian_analytic_calculated

  PUBLIC equations_matrices_create_finish,equations_matrices_create_start

  PUBLIC equationsmatrices_createfinish,equationsmatrices_createstart

  PUBLIC equations_matrices_destroy

  PUBLIC equationsmatrices_destroy

  PUBLIC equations_matrices_dynamic_lumping_type_set,equations_matrices_dynamic_storage_type_set

  PUBLIC equationsmatrices_dynamiclumpingtypeset,equationsmatrices_dynamicstoragetypeset,equationsmatrices_dynamicstructuretypeset

  !!TODO check if the elements should be create/destroy rather than initialise/finalise
  PUBLIC equations_matrices_element_initialise,equations_matrices_element_finalise

  PUBLIC equationsmatrices_elementinitialise,equationsmatrices_elementfinalise
  
  PUBLIC equations_matrices_element_add,equations_matrices_jacobian_element_add,equations_matrices_element_calculate, &
    & equations_matrices_values_initialise

  PUBLIC equationsmatrices_elementadd,equationsmatrices_jacobianelementadd,equationsmatrices_elementcalculate, &
    & equationsmatrices_valuesinitialise

  PUBLIC equations_matrices_element_matrix_finalise

  PUBLIC equationsmatrices_elementmatrixfinalise,equationsmatrices_elementmatrixinitialise
  
  PUBLIC equations_matrices_element_matrix_calculate,equations_matrices_element_matrix_setup

  PUBLIC equationsmatrices_elementmatrixcalculate,equationsmatrices_elementmatrixsetup

  PUBLIC equations_matrices_element_vector_finalise

  PUBLIC equationsmatrices_elementvectorfinalise,equationsmatrices_elementvectorinitialise

  PUBLIC equations_matrices_element_vector_calculate,equations_matrices_element_vector_setup

  PUBLIC equationsmatrices_elementvectorcalculate,equationsmatrices_elementvectorsetup

  PUBLIC equationsmatrices_jacobiantypesset

  PUBLIC equations_matrices_linear_storage_type_set

  PUBLIC equationsmatrices_linearstoragetypeset,equationsmatrices_linearstructuretypeset

  PUBLIC equationsmatrices_nodalinitialise,equationsmatrices_nodalfinalise

  PUBLIC equationsmatrices_nodeadd,equationsmatrices_jacobiannodeadd,equationsmatrices_nodalcalculate

  PUBLIC equationsmatrices_nodalmatrixfinalise,equationsmatrices_nodalmatrixinitialise

  PUBLIC equationsmatrices_nodalmatrixcalculate,equationsmatrices_nodalmatrixsetup

  PUBLIC equationsmatrices_nodalvectorfinalise,equationsmatrices_nodalvectorinitialise

  PUBLIC equationsmatrices_nodalvectorcalculate,equationsmatrices_nodalvectorsetup

  PUBLIC equationsmatrices_nonlinearstoragetypeset,equationsmatrices_nonlinearstructuretypeset

  PUBLIC equations_matrices_output,equations_matrices_jacobian_output

  PUBLIC equationsmatrices_output,equationsmatrices_jacobianoutput

CONTAINS

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

  SUBROUTINE equations_jacobian_finalise(EQUATIONS_JACOBIAN,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_jacobian_type), POINTER :: EQUATIONS_JACOBIAN
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    
    enters("EQUATIONS_JACOBIAN_FINALISE",err,error,*999)

    IF(ASSOCIATED(equations_jacobian)) THEN
      IF(ASSOCIATED(equations_jacobian%JACOBIAN)) CALL distributed_matrix_destroy(equations_jacobian%JACOBIAN,err,error,*999)
      CALL equations_matrices_element_matrix_finalise(equations_jacobian%ELEMENT_JACOBIAN,err,error,*999)
      CALL equationsmatrices_nodalmatrixfinalise(equations_jacobian%NodalJacobian,err,error,*999)
    ENDIF
    
    exits("EQUATIONS_JACOBIAN_FINALISE")
    RETURN
999 errorsexits("EQUATIONS_JACOBIAN_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_jacobian_finalise

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

  SUBROUTINE equations_jacobian_initialise(NONLINEAR_MATRICES,MATRIX_NUMBER,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
    INTEGER(INTG), INTENT(IN) :: MATRIX_NUMBER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DUMMY_ERR
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    TYPE(equations_mapping_nonlinear_type), POINTER :: NONLINEAR_MAPPING
    TYPE(varying_string) :: DUMMY_ERROR

    enters("EQUATIONS_JACOBIAN_INITIALISE",err,error,*998)
 
    IF(ASSOCIATED(nonlinear_matrices)) THEN
      equations_matrices=>nonlinear_matrices%EQUATIONS_MATRICES
      IF(ASSOCIATED(equations_matrices)) THEN
        equations_mapping=>equations_matrices%EQUATIONS_MAPPING
        IF(ASSOCIATED(equations_mapping)) THEN
          nonlinear_mapping=>equations_mapping%NONLINEAR_MAPPING
          IF(ASSOCIATED(nonlinear_mapping)) THEN
            IF(ALLOCATED(nonlinear_matrices%JACOBIANS)) THEN
              IF(ASSOCIATED(nonlinear_matrices%JACOBIANS(matrix_number)%PTR)) THEN
                CALL flagerror("Nonlinear matrices Jacobian is already associated.",err,error,*998)
              ELSE
                ALLOCATE(nonlinear_matrices%JACOBIANS(matrix_number)%PTR,stat=err)
                IF(err/=0) CALL flagerror("Could not allocate equations Jacobian.",err,error,*999)
                nonlinear_matrices%JACOBIANS(matrix_number)%PTR%JACOBIAN_NUMBER=matrix_number
                nonlinear_matrices%JACOBIANS(matrix_number)%PTR%NONLINEAR_MATRICES=>nonlinear_matrices
                nonlinear_matrices%JACOBIANS(matrix_number)%PTR%STORAGE_TYPE=matrix_block_storage_type
                nonlinear_matrices%JACOBIANS(matrix_number)%PTR%STRUCTURE_TYPE=equations_matrix_no_structure
                nonlinear_matrices%JACOBIANS(matrix_number)%PTR%NUMBER_OF_COLUMNS= &
                    & nonlinear_mapping%JACOBIAN_TO_VAR_MAP(matrix_number)%NUMBER_OF_COLUMNS
                nonlinear_matrices%JACOBIANS(matrix_number)%PTR%UPDATE_JACOBIAN=.true.
                nonlinear_matrices%JACOBIANS(matrix_number)%PTR%FIRST_ASSEMBLY=.true.
                nonlinear_mapping%JACOBIAN_TO_VAR_MAP(matrix_number)%JACOBIAN=>nonlinear_matrices%JACOBIANS(matrix_number)%PTR
                NULLIFY(nonlinear_matrices%JACOBIANS(matrix_number)%PTR%JACOBIAN)
                CALL equationsmatrices_elementmatrixinitialise(nonlinear_matrices%JACOBIANS(matrix_number)%PTR% &
                    & element_jacobian,err,error,*999)
                CALL equationsmatrices_nodalmatrixinitialise(nonlinear_matrices%JACOBIANS(matrix_number)%PTR% &
                    & nodaljacobian,err,error,*999)
                nonlinear_matrices%JACOBIANS(matrix_number)%PTR%JACOBIAN_CALCULATION_TYPE= &
                  & equations_jacobian_finite_difference_calculated
              ENDIF
            ELSE
              CALL flagerror("Equations matrices nonlinear matrieces Jacobian is not allocated.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Equations mapping nonlinear mapping is not associated.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Equations mapping is not associated.",err,error,*998)
        ENDIF
      ELSE
        CALL flagerror("Nonlinear matrices equations matrices is not associated.",err,error,*998)
      ENDIF
    ELSE
      CALL flagerror("Nonlinear matrices is not associated.",err,error,*998)
    ENDIF
    
    exits("EQUATIONS_JACOBIAN_INITIALISE")
    RETURN
999 CALL equations_jacobian_finalise(nonlinear_matrices%JACOBIANS(matrix_number)%PTR,dummy_err,dummy_error,*998)
998 errorsexits("EQUATIONS_JACOBIAN_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_jacobian_initialise

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

  SUBROUTINE equations_matrices_create_finish(EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DUMMY_ERR,matrix_idx,NUMBER_OF_NON_ZEROS
    INTEGER(INTG), POINTER :: ROW_INDICES(:),COLUMN_INDICES(:)
    TYPE(domain_mapping_type), POINTER :: ROW_DOMAIN_MAP,COLUMN_DOMAIN_MAP
    TYPE(equations_jacobian_type), POINTER :: JACOBIAN_MATRIX
    TYPE(equations_mapping_dynamic_type), POINTER :: DYNAMIC_MAPPING
    TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
    TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
    TYPE(equations_matrices_source_type), POINTER :: SOURCE_VECTOR
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    TYPE(equations_mapping_nonlinear_type), POINTER :: NONLINEAR_MAPPING
    TYPE(varying_string) :: DUMMY_ERROR,LOCAL_ERROR
    type(linkedlist),pointer :: list(:) 
    NULLIFY(row_indices)
    NULLIFY(column_indices)

    enters("EQUATIONS_MATRICES_CREATE_FINISH",err,error,*998)

    IF(ASSOCIATED(equations_matrices)) THEN
      IF(equations_matrices%EQUATIONS_MATRICES_FINISHED) THEN
        CALL flagerror("Equations matrices have already been finished.",err,error,*998)
      ELSE
        equations_mapping=>equations_matrices%EQUATIONS_MAPPING
        IF(ASSOCIATED(equations_mapping)) THEN
          row_domain_map=>equations_mapping%ROW_DOFS_MAPPING
          IF(ASSOCIATED(row_domain_map)) THEN
            dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
            IF(ASSOCIATED(dynamic_matrices)) THEN
              !Dynamic matrices
              dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
              IF(ASSOCIATED(dynamic_mapping)) THEN
                !Now create the individual dynamic equations matrices
                DO matrix_idx=1,dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES
                  equations_matrix=>dynamic_matrices%MATRICES(matrix_idx)%PTR
                  IF(ASSOCIATED(equations_matrix)) THEN
                    column_domain_map=>dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrix_idx)%COLUMN_DOFS_MAPPING
                    IF(ASSOCIATED(column_domain_map)) THEN
                      !Create the distributed equations matrix
                      CALL distributed_matrix_create_start(row_domain_map,column_domain_map,equations_matrices% &
                           & dynamic_matrices%MATRICES(matrix_idx)%PTR%MATRIX,err,error,*999)
                      CALL distributed_matrix_data_type_set(equations_matrix%MATRIX,matrix_vector_dp_type,err,error,*999)
                      CALL distributed_matrix_storage_type_set(equations_matrix%MATRIX,equations_matrix%STORAGE_TYPE,err,error,*999)
                      !Calculate and set the matrix structure/sparsity pattern
                      IF(equations_matrix%STORAGE_TYPE/=distributed_matrix_block_storage_type.AND. &
                        & equations_matrix%STORAGE_TYPE/=distributed_matrix_diagonal_storage_type) THEN
                        CALL equationsmatrix_structurecalculate(equations_matrix,number_of_non_zeros,row_indices,column_indices, &
                          & list,err,error,*999)
                        CALL distributed_matrix_linklist_set(equations_matrix%MATRIX,list,err,error,*999)
                        CALL distributed_matrix_number_non_zeros_set(equations_matrix%MATRIX,number_of_non_zeros,err,error,*999)
                        CALL distributed_matrix_storage_locations_set(equations_matrix%MATRIX,row_indices,column_indices, &
                          & err,error,*999)
                        IF(ASSOCIATED(row_indices)) DEALLOCATE(row_indices)
                        IF(ASSOCIATED(column_indices)) DEALLOCATE(column_indices)
                      ENDIF
                      CALL distributed_matrix_create_finish(equations_matrix%MATRIX,err,error,*999)
                    ELSE
                      local_error="Column domain map for dynamic matrix number "// &
                        & trim(numbertovstring(matrix_idx,"*",err,error))//" is not associated."
                      CALL flagerror(local_error,err,error,*999)
                    ENDIF
                  ELSE
                    local_error="Equations matrix for dynamic matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))// &
                      & " is not associated."
                    CALL flagerror(local_error,err,error,*999)
                  ENDIF
                ENDDO !matrix_idx                
              ELSE
                CALL flagerror("Equations mapping dynamic mapping is not associated.",err,error,*999)                
              ENDIF
            ENDIF
            linear_matrices=>equations_matrices%LINEAR_MATRICES
            IF(ASSOCIATED(linear_matrices)) THEN
              !Linear matrices
              linear_mapping=>equations_mapping%LINEAR_MAPPING
              IF(ASSOCIATED(linear_mapping)) THEN
                !Now create the individual linear equations matrices
                DO matrix_idx=1,linear_matrices%NUMBER_OF_LINEAR_MATRICES
                  equations_matrix=>linear_matrices%MATRICES(matrix_idx)%PTR
                  IF(ASSOCIATED(equations_matrix)) THEN
                    column_domain_map=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrix_idx)%COLUMN_DOFS_MAPPING
                    IF(ASSOCIATED(column_domain_map)) THEN
                      !Create the distributed equations matrix
                      CALL distributed_matrix_create_start(row_domain_map,column_domain_map,equations_matrices% &
                           & linear_matrices%MATRICES(matrix_idx)%PTR%MATRIX,err,error,*999)
                      CALL distributed_matrix_data_type_set(equations_matrix%MATRIX,matrix_vector_dp_type,err,error,*999)
                      CALL distributed_matrix_storage_type_set(equations_matrix%MATRIX,equations_matrix%STORAGE_TYPE,err,error,*999)
                      !Calculate and set the matrix structure/sparsity pattern
                      IF(equations_matrix%STORAGE_TYPE/=distributed_matrix_block_storage_type.AND. &
                        & equations_matrix%STORAGE_TYPE/=distributed_matrix_diagonal_storage_type) THEN
                        CALL equationsmatrix_structurecalculate(equations_matrix,number_of_non_zeros,row_indices,column_indices, &
                          & list,err,error,*999)
                        CALL distributed_matrix_linklist_set(equations_matrix%MATRIX,list,err,error,*999)
                        CALL distributed_matrix_number_non_zeros_set(equations_matrix%MATRIX,number_of_non_zeros,err,error,*999)
                        CALL distributed_matrix_storage_locations_set(equations_matrix%MATRIX,row_indices,column_indices, &
                          & err,error,*999)
                        IF(ASSOCIATED(row_indices)) DEALLOCATE(row_indices)
                        IF(ASSOCIATED(column_indices)) DEALLOCATE(column_indices)
                      ENDIF
                      CALL distributed_matrix_create_finish(equations_matrix%MATRIX,err,error,*999)
                    ELSE
                      local_error="Column domain map for linear matrix number "// &
                        & trim(numbertovstring(matrix_idx,"*",err,error))//" is not associated."
                      CALL flagerror(local_error,err,error,*999)
                    ENDIF
                  ELSE
                    local_error="Equations matrix for linear matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))// &
                      & " is not associated."
                    CALL flagerror(local_error,err,error,*999)
                  ENDIF
                ENDDO !matrix_idx
              ELSE
                CALL flagerror("Equations mapping linear mapping is not associated.",err,error,*999)                
              ENDIF
            ENDIF
            nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
            IF(ASSOCIATED(nonlinear_matrices)) THEN
              !Nonlinear matrices
              nonlinear_mapping=>equations_mapping%NONLINEAR_MAPPING
              IF(ASSOCIATED(nonlinear_mapping)) THEN
                !Set up the Jacobian matrices
                DO matrix_idx=1,nonlinear_matrices%NUMBER_OF_JACOBIANS
                  jacobian_matrix=>nonlinear_matrices%JACOBIANS(matrix_idx)%PTR
                  IF(ASSOCIATED(jacobian_matrix)) THEN
                    column_domain_map=>nonlinear_mapping%JACOBIAN_TO_VAR_MAP(matrix_idx)%COLUMN_DOFS_MAPPING
                    IF(ASSOCIATED(column_domain_map)) THEN
  !!TODO: Set the distributed matrix not to allocate the data if the Jacobian is not calculated.
                      !Create the distributed Jacobian matrix
                      CALL distributed_matrix_create_start(row_domain_map,column_domain_map,jacobian_matrix%JACOBIAN,err,error,*999)
                      CALL distributed_matrix_data_type_set(jacobian_matrix%JACOBIAN,matrix_vector_dp_type,err,error,*999)
                      CALL distributed_matrix_storage_type_set(jacobian_matrix%JACOBIAN,jacobian_matrix%STORAGE_TYPE,err,error,*999)
                      !Calculate and set the matrix structure/sparsity pattern
                      IF(jacobian_matrix%STORAGE_TYPE/=distributed_matrix_block_storage_type.AND. &
                        & jacobian_matrix%STORAGE_TYPE/=distributed_matrix_diagonal_storage_type) THEN
                        CALL jacobianmatrix_structurecalculate(jacobian_matrix,number_of_non_zeros,row_indices,column_indices, &
                          & err,error,*999)
                        CALL distributed_matrix_number_non_zeros_set(jacobian_matrix%JACOBIAN,number_of_non_zeros,err,error,*999)
                        CALL distributed_matrix_storage_locations_set(jacobian_matrix%JACOBIAN,row_indices,column_indices, &
                          & err,error,*999)
                        IF(ASSOCIATED(row_indices)) DEALLOCATE(row_indices)
                        IF(ASSOCIATED(column_indices)) DEALLOCATE(column_indices)
                      ENDIF
                      CALL distributed_matrix_create_finish(jacobian_matrix%JACOBIAN,err,error,*999)
                    ELSE
                      CALL flagerror("Column domain map is not associated.",err,error,*999)
                    ENDIF
                  ELSE
                    local_error="Jacobian matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))//" is not associated."
                    CALL flagerror(local_error,err,error,*999)
                  ENDIF
                ENDDO
                !Set up the residual vector                
                CALL distributed_vector_create_start(row_domain_map,equations_matrices%NONLINEAR_MATRICES%RESIDUAL,err,error,*999)
                CALL distributed_vector_data_type_set(nonlinear_matrices%RESIDUAL,matrix_vector_dp_type,err,error,*999)
                CALL distributed_vector_create_finish(nonlinear_matrices%RESIDUAL,err,error,*999)
                !Initialise the residual vector to zero for time dependent problems so that the previous residual is set to zero
                CALL distributed_vector_all_values_set(nonlinear_matrices%RESIDUAL,0.0_dp,err,error,*999)
              ELSE
                CALL flagerror("Equations mapping nonlinear mapping is not associated.",err,error,*999)
              ENDIF
            ENDIF
            rhs_vector=>equations_matrices%RHS_VECTOR
            IF(ASSOCIATED(rhs_vector)) THEN
              !Set up the equations RHS vector          
              CALL distributed_vector_create_start(row_domain_map,equations_matrices%RHS_VECTOR%VECTOR,err,error,*999)
              CALL distributed_vector_data_type_set(rhs_vector%VECTOR,matrix_vector_dp_type,err,error,*999)
              CALL distributed_vector_create_finish(rhs_vector%VECTOR,err,error,*999)
            ENDIF
            source_vector=>equations_matrices%SOURCE_VECTOR
            IF(ASSOCIATED(source_vector)) THEN
              !Set up the equations source vector          
              CALL distributed_vector_create_start(row_domain_map,equations_matrices%SOURCE_VECTOR%VECTOR,err,error,*999)
              CALL distributed_vector_data_type_set(source_vector%VECTOR,matrix_vector_dp_type,err,error,*999)
              CALL distributed_vector_create_finish(source_vector%VECTOR,err,error,*999)
            ENDIF
            !Finish up
            equations_matrices%EQUATIONS_MATRICES_FINISHED=.true.
          ELSE
            CALL flagerror("Row domain map is not associated.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Equations mapping is not associated.",err,error,*998)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*998)
    ENDIF
       
    exits("EQUATIONS_MATRICES_CREATE_FINISH")
    RETURN
999 IF(ASSOCIATED(row_indices)) DEALLOCATE(row_indices)
    IF(ASSOCIATED(column_indices)) DEALLOCATE(column_indices)
    CALL equations_matrices_finalise(equations_matrices,dummy_err,dummy_error,*998)
998 errorsexits("EQUATIONS_MATRICES_CREATE_FINISH",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_create_finish

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

  SUBROUTINE equations_matrices_create_start(EQUATIONS,EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_type), POINTER :: EQUATIONS
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DUMMY_ERR
    TYPE(varying_string) :: DUMMY_ERROR    

    enters("EQUATIONS_MATRICES_CREATE_START",err,error,*998)

    IF(ASSOCIATED(equations)) THEN      
      IF(equations%EQUATIONS_FINISHED) THEN
        IF(ASSOCIATED(equations_matrices)) THEN
          CALL flagerror("Equations matrices is already associated.",err,error,*998)
        ELSE
          NULLIFY(equations_matrices)
          !Initialise the equations matrices
          CALL equations_matrices_initialise(equations,err,error,*999)
          equations_matrices=>equations%EQUATIONS_MATRICES
        ENDIF
      ELSE
        CALL flagerror("Equations has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations is not associated.",err,error,*998)
    ENDIF
    
    exits("EQUATIONS_MATRICES_CREATE_START")
    RETURN
999 CALL equations_matrices_finalise(equations%EQUATIONS_MATRICES,dummy_err,dummy_error,*998)
998 errorsexits("EQUATIONS_MATRICES_CREATE_START",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_create_start

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

  SUBROUTINE equations_matrices_destroy(EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(equations_matrices)) THEN
      CALL equations_matrices_finalise(equations_matrices,err,error,*999)
    ELSE
      CALL flagerror("Equations matrices is not associated",err,error,*999)
    ENDIF
        
    exits("EQUATIONS_MATRICES_DESTROY")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_DESTROY",err,error)    
    RETURN 1
   
  END SUBROUTINE equations_matrices_destroy

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

  SUBROUTINE equations_matrices_element_matrix_calculate(ELEMENT_MATRIX,UPDATE_MATRIX,ROW_ELEMENT_NUMBERS,COLUMN_ELEMENT_NUMBERS, &
    & rows_field_variable,cols_field_variable,err,error,*)

    !Argument variables
    TYPE(element_matrix_type) :: ELEMENT_MATRIX
    LOGICAL :: UPDATE_MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_ELEMENT_NUMBERS(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_ELEMENT_NUMBERS(:)
    TYPE(field_variable_type), POINTER :: ROWS_FIELD_VARIABLE
    TYPE(field_variable_type), POINTER :: COLS_FIELD_VARIABLE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: component_idx,derivative,derivative_idx,global_ny,local_ny,node,node_idx,version,dataPointIdx, &
      & localDataPointNumber,elementIdx,rowElementNumber,colElementNumber
    TYPE(basis_type), POINTER :: BASIS
    TYPE(domain_elements_type), POINTER :: ELEMENTS_TOPOLOGY
    TYPE(decompositiondatapointstype), POINTER :: decompositionData
    TYPE(varying_string) :: LOCAL_ERROR

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

    IF(ASSOCIATED(rows_field_variable)) THEN
      IF(ASSOCIATED(cols_field_variable)) THEN
        element_matrix%NUMBER_OF_ROWS=0
        element_matrix%NUMBER_OF_COLUMNS=0
        IF(update_matrix) THEN
          IF(ASSOCIATED(rows_field_variable,cols_field_variable)) THEN
            !Row and columns variable is the same.
            DO component_idx=1,rows_field_variable%NUMBER_OF_COMPONENTS
              elements_topology=>rows_field_variable%COMPONENTS(component_idx)%DOMAIN%TOPOLOGY%ELEMENTS
              DO elementidx=1,SIZE(row_element_numbers)
                rowelementnumber=row_element_numbers(elementidx)
                IF(rowelementnumber>=1.AND.rowelementnumber<=elements_topology%TOTAL_NUMBER_OF_ELEMENTS) THEN
                  SELECT CASE(rows_field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE)
                  CASE(field_constant_interpolation)
                    local_ny=rows_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
                    global_ny=rows_field_variable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(local_ny)
                    element_matrix%NUMBER_OF_ROWS=element_matrix%NUMBER_OF_ROWS+1
                    element_matrix%NUMBER_OF_COLUMNS=element_matrix%NUMBER_OF_COLUMNS+1
                    element_matrix%ROW_DOFS(element_matrix%NUMBER_OF_ROWS)=local_ny
                    element_matrix%COLUMN_DOFS(element_matrix%NUMBER_OF_COLUMNS)=global_ny
                  CASE(field_element_based_interpolation)
                    local_ny=rows_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP% &
                      & elements(rowelementnumber)
                    global_ny=rows_field_variable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(local_ny)
                    element_matrix%NUMBER_OF_ROWS=element_matrix%NUMBER_OF_ROWS+1
                    element_matrix%NUMBER_OF_COLUMNS=element_matrix%NUMBER_OF_COLUMNS+1
                    element_matrix%ROW_DOFS(element_matrix%NUMBER_OF_ROWS)=local_ny
                    element_matrix%COLUMN_DOFS(element_matrix%NUMBER_OF_COLUMNS)=global_ny
                  CASE(field_node_based_interpolation)
                    basis=>elements_topology%ELEMENTS(rowelementnumber)%BASIS
                    DO node_idx=1,basis%NUMBER_OF_NODES
                      node=elements_topology%ELEMENTS(rowelementnumber)%ELEMENT_NODES(node_idx)
                      DO derivative_idx=1,basis%NUMBER_OF_DERIVATIVES(node_idx)
                        derivative=elements_topology%ELEMENTS(rowelementnumber)%ELEMENT_DERIVATIVES(derivative_idx,node_idx)
                        version=elements_topology%ELEMENTS(rowelementnumber)%elementVersions(derivative_idx,node_idx)
                        local_ny=rows_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node)% &
                          & derivatives(derivative)%VERSIONS(version)
                        global_ny=rows_field_variable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(local_ny)
                        element_matrix%NUMBER_OF_ROWS=element_matrix%NUMBER_OF_ROWS+1
                        element_matrix%NUMBER_OF_COLUMNS=element_matrix%NUMBER_OF_COLUMNS+1
                        element_matrix%ROW_DOFS(element_matrix%NUMBER_OF_ROWS)=local_ny
                        element_matrix%COLUMN_DOFS(element_matrix%NUMBER_OF_COLUMNS)=global_ny
                      ENDDO !derivative_idx
                    ENDDO !node_idx
                  CASE(field_grid_point_based_interpolation)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(field_gauss_point_based_interpolation)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(field_data_point_based_interpolation)
                    decompositiondata=>rows_field_variable%COMPONENTS(component_idx)%DOMAIN%DECOMPOSITION%TOPOLOGY%dataPoints
                    DO datapointidx=1,decompositiondata%elementDataPoint(rowelementnumber)%numberOfProjectedData
                      localdatapointnumber=decompositiondata%elementDataPoint(rowelementnumber)% &
                        & dataindices(datapointidx)%localNumber
                      local_ny=rows_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%DATA_POINT_PARAM2DOF_MAP% &
                        & data_points(localdatapointnumber)
                      global_ny=rows_field_variable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(local_ny)
                      element_matrix%NUMBER_OF_ROWS=element_matrix%NUMBER_OF_ROWS+1
                      element_matrix%NUMBER_OF_COLUMNS=element_matrix%NUMBER_OF_COLUMNS+1
                      element_matrix%ROW_DOFS(element_matrix%NUMBER_OF_ROWS)=local_ny
                      element_matrix%COLUMN_DOFS(element_matrix%NUMBER_OF_COLUMNS)=global_ny
                    ENDDO
                  CASE DEFAULT
                    local_error="The interpolation type of "// &
                      & trim(numbertovstring(rows_field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE,"*",err,error))// &
                      & " is invalid for component number "// &
                      & trim(numbertovstring(component_idx,"*",err,error))// &
                      & " of rows field variable type "// &
                      & trim(numbertovstring(rows_field_variable%VARIABLE_TYPE,"*",err,error))//"."
                    CALL flagerror(local_error,err,error,*999)          
                  END SELECT
                ELSE
                  local_error="Element number "//trim(numbertovstring(rowelementnumber,"*",err,error))// &
                    & " is invalid for component number "//trim(numbertovstring(component_idx,"*",err,error))// &
                    & " of rows field variable type "// &
                    & trim(numbertovstring(rows_field_variable%VARIABLE_TYPE,"*",err,error))// &
                    & ". The element number must be between 1 and "// &
                    & trim(numbertovstring(elements_topology%TOTAL_NUMBER_OF_ELEMENTS,"*",err,error))//"."
                  CALL flagerror(local_error,err,error,*999)
                ENDIF
              ENDDO !elementIdx
            ENDDO !component_idx
          ELSE
            !Row and column variables are different
            !Row mapping
            DO component_idx=1,rows_field_variable%NUMBER_OF_COMPONENTS
              elements_topology=>rows_field_variable%COMPONENTS(component_idx)%DOMAIN%TOPOLOGY%ELEMENTS
              DO elementidx=1,SIZE(row_element_numbers)
                rowelementnumber=row_element_numbers(elementidx)
                IF(rowelementnumber>=1.AND.rowelementnumber<=elements_topology%TOTAL_NUMBER_OF_ELEMENTS) THEN
                  SELECT CASE(rows_field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE)
                  CASE(field_constant_interpolation)
                    local_ny=rows_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
                    element_matrix%NUMBER_OF_ROWS=element_matrix%NUMBER_OF_ROWS+1
                    element_matrix%ROW_DOFS(element_matrix%NUMBER_OF_ROWS)=local_ny
                  CASE(field_element_based_interpolation)
                    local_ny=rows_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP% &
                      & elements(rowelementnumber)
                    element_matrix%NUMBER_OF_ROWS=element_matrix%NUMBER_OF_ROWS+1
                    element_matrix%ROW_DOFS(element_matrix%NUMBER_OF_ROWS)=local_ny
                  CASE(field_node_based_interpolation)
                    basis=>elements_topology%ELEMENTS(rowelementnumber)%BASIS
                    DO node_idx=1,basis%NUMBER_OF_NODES
                      node=elements_topology%ELEMENTS(rowelementnumber)%ELEMENT_NODES(node_idx)
                      DO derivative_idx=1,basis%NUMBER_OF_DERIVATIVES(node_idx)
                        derivative=elements_topology%ELEMENTS(rowelementnumber)%ELEMENT_DERIVATIVES(derivative_idx,node_idx)
                        version=elements_topology%ELEMENTS(rowelementnumber)%elementVersions(derivative_idx,node_idx)
                        local_ny=rows_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node)% &
                          & derivatives(derivative)%VERSIONS(version)
                        element_matrix%NUMBER_OF_ROWS=element_matrix%NUMBER_OF_ROWS+1
                        element_matrix%ROW_DOFS(element_matrix%NUMBER_OF_ROWS)=local_ny
                      ENDDO !derivative_idx
                    ENDDO !node_idx
                  CASE(field_grid_point_based_interpolation)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(field_gauss_point_based_interpolation)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(field_data_point_based_interpolation)
                    decompositiondata=>rows_field_variable%COMPONENTS(component_idx)%DOMAIN%DECOMPOSITION%TOPOLOGY%dataPoints
                    DO datapointidx=1,decompositiondata%elementDataPoint(colelementnumber)%numberOfProjectedData
                      localdatapointnumber=decompositiondata%elementDataPoint(colelementnumber)% &
                        & dataindices(datapointidx)%localNumber
                      local_ny=rows_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%DATA_POINT_PARAM2DOF_MAP% &
                        & data_points(localdatapointnumber)
                      global_ny=rows_field_variable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(local_ny)
                      element_matrix%NUMBER_OF_COLUMNS=element_matrix%NUMBER_OF_COLUMNS+1
                      element_matrix%COLUMN_DOFS(element_matrix%NUMBER_OF_COLUMNS)=global_ny
                    ENDDO
                  CASE DEFAULT
                    local_error="The interpolation type of "// &
                      & trim(numbertovstring(rows_field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE,"*",err,error))// &
                      & " is invalid for component number "// &
                      & trim(numbertovstring(component_idx,"*",err,error))// &
                      & " of rows field variable type "// &
                      & trim(numbertovstring(rows_field_variable%VARIABLE_TYPE,"*",err,error))//"."
                    CALL flagerror(local_error,err,error,*999)          
                  END SELECT
                ELSE
                  local_error="Row element number "//trim(numbertovstring(rowelementnumber,"*",err,error))// &
                    & " is invalid for component number "//trim(numbertovstring(component_idx,"*",err,error))// &
                    & " of rows field variable type "// &
                    & trim(numbertovstring(rows_field_variable%VARIABLE_TYPE,"*",err,error))// &
                    & ". The element number must be between 1 and "// &
                    & trim(numbertovstring(elements_topology%TOTAL_NUMBER_OF_ELEMENTS,"*",err,error))//"."
                  CALL flagerror(local_error,err,error,*999)
                ENDIF
              ENDDO !elementIdx
            ENDDO !component_idx
            !Column mapping
            DO component_idx=1,cols_field_variable%NUMBER_OF_COMPONENTS
              elements_topology=>cols_field_variable%COMPONENTS(component_idx)%DOMAIN%TOPOLOGY%ELEMENTS
              DO elementidx=1,SIZE(column_element_numbers)
                colelementnumber=column_element_numbers(elementidx)
                IF(colelementnumber>=1.AND.colelementnumber<=elements_topology%TOTAL_NUMBER_OF_ELEMENTS) THEN
                  SELECT CASE(cols_field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE)
                  CASE(field_constant_interpolation)
                    local_ny=cols_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
                    global_ny=cols_field_variable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(local_ny)
                    element_matrix%NUMBER_OF_COLUMNS=element_matrix%NUMBER_OF_COLUMNS+1
                    element_matrix%COLUMN_DOFS(element_matrix%NUMBER_OF_COLUMNS)=global_ny
                  CASE(field_element_based_interpolation)
                    local_ny=cols_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP% &
                      & elements(colelementnumber)
                    global_ny=cols_field_variable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(local_ny)
                    element_matrix%NUMBER_OF_COLUMNS=element_matrix%NUMBER_OF_COLUMNS+1
                    element_matrix%COLUMN_DOFS(element_matrix%NUMBER_OF_COLUMNS)=global_ny
                  CASE(field_node_based_interpolation)
                    basis=>elements_topology%ELEMENTS(colelementnumber)%BASIS
                    DO node_idx=1,basis%NUMBER_OF_NODES
                      node=elements_topology%ELEMENTS(colelementnumber)%ELEMENT_NODES(node_idx)
                      DO derivative_idx=1,basis%NUMBER_OF_DERIVATIVES(node_idx)
                        derivative=elements_topology%ELEMENTS(colelementnumber)%ELEMENT_DERIVATIVES(derivative_idx,node_idx)
                        version=elements_topology%ELEMENTS(colelementnumber)%elementVersions(derivative_idx,node_idx)
                        local_ny=cols_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node)% &
                          & derivatives(derivative)%VERSIONS(version)
                        global_ny=cols_field_variable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(local_ny)
                        element_matrix%NUMBER_OF_COLUMNS=element_matrix%NUMBER_OF_COLUMNS+1
                        element_matrix%COLUMN_DOFS(element_matrix%NUMBER_OF_COLUMNS)=global_ny
                      ENDDO !derivative_idx
                    ENDDO !node_idx
                  CASE(field_grid_point_based_interpolation)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(field_gauss_point_based_interpolation)
                    CALL flagerror("Not implemented.",err,error,*999)
                  CASE(field_data_point_based_interpolation)
                    decompositiondata=>cols_field_variable%COMPONENTS(component_idx)%DOMAIN%DECOMPOSITION%TOPOLOGY%dataPoints
                    DO datapointidx=1,decompositiondata%elementDataPoint(colelementnumber)%numberOfProjectedData
                      localdatapointnumber=decompositiondata%elementDataPoint(colelementnumber)% &
                        & dataindices(datapointidx)%localNumber
                      local_ny=cols_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%DATA_POINT_PARAM2DOF_MAP% &
                        & data_points(localdatapointnumber)
                      global_ny=cols_field_variable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(local_ny)
                      element_matrix%NUMBER_OF_COLUMNS=element_matrix%NUMBER_OF_COLUMNS+1
                      element_matrix%COLUMN_DOFS(element_matrix%NUMBER_OF_COLUMNS)=global_ny
                    ENDDO
                  CASE DEFAULT
                    local_error="The interpolation type of "// &
                      & trim(numbertovstring(cols_field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE,"*",err,error))// &
                      & " is invalid for component number "// &
                      & trim(numbertovstring(component_idx,"*",err,error))// &
                      & " of column field variable type "// &
                      & trim(numbertovstring(cols_field_variable%VARIABLE_TYPE,"*",err,error))//"."
                    CALL flagerror(local_error,err,error,*999)          
                  END SELECT
                ELSE
                  local_error="Column element number "//trim(numbertovstring(colelementnumber,"*",err,error))// &
                    & " is invalid for component number "//trim(numbertovstring(component_idx,"*",err,error))// &
                    & " of column field variable type "// &
                    & trim(numbertovstring(cols_field_variable%VARIABLE_TYPE,"*",err,error))// &
                    & ". The element number must be between 1 and "// &
                    & trim(numbertovstring(elements_topology%TOTAL_NUMBER_OF_ELEMENTS,"*",err,error))//"."
                  CALL flagerror(local_error,err,error,*999)
                ENDIF
              ENDDO !elementIdx
            ENDDO !component_idx
          ENDIF
          element_matrix%MATRIX=0.0_dp
        ENDIF
      ELSE
        CALL flagerror("Columns field variable is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Rows field variable is not associated.",err,error,*999)
    ENDIF
    
    exits("EQUATIONS_MATRICES_ELEMENT_MATRIX_CALCULATE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_ELEMENT_MATRIX_CALCULATE",err,error)
    RETURN 1
    
  END SUBROUTINE equations_matrices_element_matrix_calculate

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

  SUBROUTINE equations_matrices_element_matrix_finalise(ELEMENT_MATRIX,ERR,ERROR,*)

    !Argument variables
    TYPE(element_matrix_type):: ELEMENT_MATRIX
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    
    enters("EQUATIONS_MATRICES_ELEMENT_MATRIX_FINALISE",err,error,*999)

    element_matrix%MAX_NUMBER_OF_ROWS=0
    element_matrix%MAX_NUMBER_OF_COLUMNS=0
    IF(ALLOCATED(element_matrix%ROW_DOFS)) DEALLOCATE(element_matrix%ROW_DOFS)
    IF(ALLOCATED(element_matrix%COLUMN_DOFS)) DEALLOCATE(element_matrix%COLUMN_DOFS)
    IF(ALLOCATED(element_matrix%MATRIX)) DEALLOCATE(element_matrix%MATRIX)
    
    exits("EQUATIONS_MATRICES_ELEMENT_MATRIX_FINALISE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_ELEMENT_MATRIX_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_element_matrix_finalise

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

  SUBROUTINE equationsmatrices_elementmatrixinitialise(ELEMENT_MATRIX,ERR,ERROR,*)

    !Argument variables
    TYPE(element_matrix_type) :: ELEMENT_MATRIX !The element matrix to initialise
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    element_matrix%EQUATIONS_MATRIX_NUMBER=0
    element_matrix%NUMBER_OF_ROWS=0
    element_matrix%NUMBER_OF_COLUMNS=0
    element_matrix%MAX_NUMBER_OF_ROWS=0
    element_matrix%MAX_NUMBER_OF_COLUMNS=0
       
    exits("EquationsMatrices_ElementMatrixInitialise")
    RETURN
999 errorsexits("EquationsMatrices_ElementMatrixInitialise",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_elementmatrixinitialise

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

  SUBROUTINE equations_matrices_element_matrix_setup(elementMatrix,rowsFieldVariable,columnsFieldVariable, &
    & rowsnumberofelements,colsnumberofelements,err,error,*)

    !Argument variables
    TYPE(element_matrix_type) :: elementMatrix
    TYPE(field_variable_type), POINTER :: rowsFieldVariable
    TYPE(field_variable_type), POINTER :: columnsFieldVariable
    INTEGER(INTG), INTENT(IN)  :: rowsNumberOfElements
    INTEGER(INTG), INTENT(IN)  :: colsNumberOfElements
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: dummyErr, componentIdx
    TYPE(varying_string) :: dummyError

    enters("EQUATIONS_MATRICES_ELEMENT_MATRIX_SETUP",err,error,*998)

    IF(ASSOCIATED(rowsfieldvariable)) THEN
      IF(ASSOCIATED(columnsfieldvariable)) THEN
        elementmatrix%MAX_NUMBER_OF_ROWS = 0
        DO componentidx=1,rowsfieldvariable%NUMBER_OF_COMPONENTS
          elementmatrix%MAX_NUMBER_OF_ROWS=elementmatrix%MAX_NUMBER_OF_ROWS+ &
            & rowsfieldvariable%COMPONENTS(componentidx)%maxNumberElementInterpolationParameters
        ENDDO
        elementmatrix%MAX_NUMBER_OF_ROWS=elementmatrix%MAX_NUMBER_OF_ROWS*rowsnumberofelements
        elementmatrix%MAX_NUMBER_OF_COLUMNS = 0
        DO componentidx=1,columnsfieldvariable%NUMBER_OF_COMPONENTS
          elementmatrix%MAX_NUMBER_OF_COLUMNS=elementmatrix%MAX_NUMBER_OF_COLUMNS+ &
            & columnsfieldvariable%COMPONENTS(componentidx)%maxNumberElementInterpolationParameters
        ENDDO
        elementmatrix%MAX_NUMBER_OF_COLUMNS=elementmatrix%MAX_NUMBER_OF_COLUMNS*colsnumberofelements
        IF(ALLOCATED(elementmatrix%ROW_DOFS)) THEN
          CALL flagerror("Element matrix row dofs already allocated.",err,error,*999)
        ELSE
          ALLOCATE(elementmatrix%ROW_DOFS(elementmatrix%MAX_NUMBER_OF_ROWS),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate element matrix row dofs.",err,error,*999)
        ENDIF
        IF(ALLOCATED(elementmatrix%COLUMN_DOFS)) THEN
          CALL flagerror("Element matrix column dofs already allocated.",err,error,*999)
        ELSE
          ALLOCATE(elementmatrix%COLUMN_DOFS(elementmatrix%MAX_NUMBER_OF_COLUMNS),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate element matrix column dofs.",err,error,*999)
        ENDIF
        IF(ALLOCATED(elementmatrix%MATRIX)) THEN
          CALL flagerror("Element matrix already allocated.",err,error,*999)
        ELSE
          ALLOCATE(elementmatrix%MATRIX(elementmatrix%MAX_NUMBER_OF_ROWS,elementmatrix%MAX_NUMBER_OF_COLUMNS),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate element matrix.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Columns field variable is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Rows field variable is not associated.",err,error,*999)
    ENDIF
    
    exits("EQUATIONS_MATRICES_ELEMENT_MATRIX_SETUP")
    RETURN
999 CALL equations_matrices_element_matrix_finalise(elementmatrix,dummyerr,dummyerror,*998)
998 errorsexits("EQUATIONS_MATRICES_ELEMENT_MATRIX_SETUP",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_element_matrix_setup

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

  SUBROUTINE equations_matrices_element_vector_calculate(ELEMENT_VECTOR,UPDATE_VECTOR,ELEMENT_NUMBER,ROWS_FIELD_VARIABLE, &
    & err,error,*)

    !Argument variables
    TYPE(element_vector_type) :: ELEMENT_VECTOR
    LOGICAL :: UPDATE_VECTOR
    INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
    TYPE(field_variable_type), POINTER :: ROWS_FIELD_VARIABLE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: component_idx,derivative,derivative_idx,local_ny,node,node_idx,version,dataPointIdx,localDataPointNumber
    TYPE(basis_type), POINTER :: BASIS
    TYPE(domain_elements_type), POINTER :: ELEMENTS_TOPOLOGY
    TYPE(decompositiondatapointstype), POINTER :: decompositionData
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("EQUATIONS_MATRICES_ELEMENT_VECTOR_CALCULATE",err,error,*999)

    IF(ASSOCIATED(rows_field_variable)) THEN
      !Calculate the rows for the element vector
      element_vector%NUMBER_OF_ROWS=0
      IF(update_vector) THEN
        DO component_idx=1,rows_field_variable%NUMBER_OF_COMPONENTS
          elements_topology=>rows_field_variable%COMPONENTS(component_idx)%DOMAIN%TOPOLOGY%ELEMENTS
          IF(element_number>=1.AND.element_number<=elements_topology%TOTAL_NUMBER_OF_ELEMENTS) THEN
            SELECT CASE(rows_field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE)
            CASE(field_constant_interpolation)
              local_ny=rows_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
              element_vector%NUMBER_OF_ROWS=element_vector%NUMBER_OF_ROWS+1
              element_vector%ROW_DOFS(element_vector%NUMBER_OF_ROWS)=local_ny
            CASE(field_element_based_interpolation)
              local_ny=rows_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP% &
                & elements(element_number)
              element_vector%NUMBER_OF_ROWS=element_vector%NUMBER_OF_ROWS+1
              element_vector%ROW_DOFS(element_vector%NUMBER_OF_ROWS)=local_ny
            CASE(field_node_based_interpolation)
              basis=>elements_topology%ELEMENTS(element_number)%BASIS
              DO node_idx=1,basis%NUMBER_OF_NODES
                node=elements_topology%ELEMENTS(element_number)%ELEMENT_NODES(node_idx)
                DO derivative_idx=1,basis%NUMBER_OF_DERIVATIVES(node_idx)
                  derivative=elements_topology%ELEMENTS(element_number)%ELEMENT_DERIVATIVES(derivative_idx,node_idx)
                  version=elements_topology%ELEMENTS(element_number)%elementVersions(derivative_idx,node_idx)
                  local_ny=rows_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node)% &
                    & derivatives(derivative)%VERSIONS(version)
                  element_vector%NUMBER_OF_ROWS=element_vector%NUMBER_OF_ROWS+1
                  element_vector%ROW_DOFS(element_vector%NUMBER_OF_ROWS)=local_ny
                ENDDO !derivative_idx
              ENDDO !node_idx
            CASE(field_grid_point_based_interpolation)
              CALL flagerror("Not implemented.",err,error,*999)
            CASE(field_gauss_point_based_interpolation)
              CALL flagerror("Not implemented.",err,error,*999)
            CASE(field_data_point_based_interpolation)
              decompositiondata=>rows_field_variable%COMPONENTS(component_idx)%DOMAIN%DECOMPOSITION%TOPOLOGY%dataPoints
              DO datapointidx=1,decompositiondata%elementDataPoint(element_number)%numberOfProjectedData
                localdatapointnumber=decompositiondata%elementDataPoint(element_number)% &
                  & dataindices(datapointidx)%localNumber
                local_ny=rows_field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%DATA_POINT_PARAM2DOF_MAP% &
                  & data_points(localdatapointnumber)
                element_vector%NUMBER_OF_ROWS=element_vector%NUMBER_OF_ROWS+1
                element_vector%ROW_DOFS(element_vector%NUMBER_OF_ROWS)=local_ny
              ENDDO
            CASE DEFAULT
              local_error="The interpolation type of "// &
                & trim(numbertovstring(rows_field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE,"*",err,error))// &
                & " is invalid for component number "// &
                & trim(numbertovstring(component_idx,"*",err,error))// &
                & " of rows field variable type "// &
                & trim(numbertovstring(rows_field_variable%VARIABLE_TYPE,"*",err,error))//"."
              CALL flagerror(local_error,err,error,*999)          
            END SELECT
          ELSE
            local_error="Element number "//trim(numbertovstring(element_number,"*",err,error))// &
              & " is invalid for component number "//trim(numbertovstring(component_idx,"*",err,error))// &
              & " of rows field variable type "//trim(numbertovstring(rows_field_variable%VARIABLE_TYPE,"*",err,error))// &
              & ". The element number must be between 1 and "// &
              & trim(numbertovstring(elements_topology%TOTAL_NUMBER_OF_ELEMENTS,"*",err,error))//"."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ENDDO !component_idx
        element_vector%VECTOR=0.0_dp
      ENDIF
    ELSE
      CALL flagerror("Rows field variable is not associated.",err,error,*999)
    ENDIF
    
    exits("EQUATIONS_MATRICES_ELEMENT_VECTOR_CALCULATE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_ELEMENT_VECTOR_CALCULATE",err,error)
    RETURN 1
    
  END SUBROUTINE equations_matrices_element_vector_calculate

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

  SUBROUTINE equations_matrices_element_vector_finalise(ELEMENT_VECTOR,ERR,ERROR,*)

    !Argument variables
    TYPE(element_vector_type):: ELEMENT_VECTOR
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    
    enters("EQUATIONS_MATRICES_ELEMENT_VECTOR_FINALISE",err,error,*999)

    IF(ALLOCATED(element_vector%ROW_DOFS)) DEALLOCATE(element_vector%ROW_DOFS)
    IF(ALLOCATED(element_vector%VECTOR)) DEALLOCATE(element_vector%VECTOR)
    
    exits("EQUATIONS_MATRICES_ELEMENT_VECTOR_FINALISE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_ELEMENT_VECTOR_FINALISE",err,error)

    RETURN 1
  END SUBROUTINE equations_matrices_element_vector_finalise

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

  SUBROUTINE equationsmatrices_elementvectorinitialise(ELEMENT_VECTOR,ERR,ERROR,*)

    !Argument variables
    TYPE(element_vector_type) :: ELEMENT_VECTOR !The element vector to initialise
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    element_vector%NUMBER_OF_ROWS=0
    element_vector%MAX_NUMBER_OF_ROWS=0
       
    exits("EquationsMatrices_ElementVectorInitialise")
    RETURN
999 errorsexits("EquationsMatrices_ElementVectorInitialise",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_elementvectorinitialise

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

  SUBROUTINE equations_matrices_element_vector_setup(elementVector,rowsFieldVariable,err,error,*)

    !Argument variables
    TYPE(element_vector_type) :: elementVector
    TYPE(field_variable_type), POINTER :: rowsFieldVariable
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: DUMMY_ERR,componentIdx
    TYPE(varying_string) :: dummyError

    enters("EQUATIONS_MATRICES_ELEMENT_VECTOR_SETUP",err,error,*998)

    IF(ASSOCIATED(rowsfieldvariable)) THEN
      elementvector%MAX_NUMBER_OF_ROWS = 0
      DO componentidx=1,rowsfieldvariable%NUMBER_OF_COMPONENTS
        elementvector%MAX_NUMBER_OF_ROWS=elementvector%MAX_NUMBER_OF_ROWS+ &
          & rowsfieldvariable%COMPONENTS(componentidx)%maxNumberElementInterpolationParameters
      ENDDO
      IF(ALLOCATED(elementvector%ROW_DOFS)) THEN
        CALL flagerror("Element vector row dofs is already allocated.",err,error,*999)
      ELSE
        ALLOCATE(elementvector%ROW_DOFS(elementvector%MAX_NUMBER_OF_ROWS),stat=err)
        IF(err/=0) CALL flagerror("Could not allocate element vector row dofs.",err,error,*999)
      ENDIF
      IF(ALLOCATED(elementvector%VECTOR)) THEN
        CALL flagerror("Element vector vector already allocated.",err,error,*999)
      ELSE
        ALLOCATE(elementvector%VECTOR(elementvector%MAX_NUMBER_OF_ROWS),stat=err)
        IF(err/=0) CALL flagerror("Could not allocate element vector vector.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Rows field variable is not associated.",err,error,*999)
    ENDIF
    
    exits("EQUATIONS_MATRICES_ELEMENT_VECTOR_SETUP")
    RETURN
999 CALL equations_matrices_element_vector_finalise(elementvector,dummy_err,dummyerror,*998)
998 errorsexits("EQUATIONS_MATRICES_ELEMENT_VECTOR_SETUP",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_element_vector_setup

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

  SUBROUTINE equations_matrices_element_add(EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: column_idx,matrix_idx,row_idx
    REAL(DP) :: SUM
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
    TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
    TYPE(equations_matrices_source_type), POINTER :: SOURCE_VECTOR
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    TYPE(varying_string) :: LOCAL_ERROR
    
#ifdef TAUPROF
    CALL tau_static_phase_start("EQUATIONS_MATRICES_ELEMENT_ADD()")
#endif

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

    IF(ASSOCIATED(equations_matrices)) THEN
      dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
      IF(ASSOCIATED(dynamic_matrices)) THEN
        !Add the element matrices
        DO matrix_idx=1,dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES
          equations_matrix=>dynamic_matrices%MATRICES(matrix_idx)%PTR
          IF(ASSOCIATED(equations_matrix)) THEN
            IF(equations_matrix%UPDATE_MATRIX) THEN
              !Handle lumped matrices
              IF(equations_matrix%LUMPED) THEN
                DO row_idx=1,equations_matrix%ELEMENT_MATRIX%NUMBER_OF_ROWS
                  sum=0.0_dp
                  DO column_idx=1,equations_matrix%ELEMENT_MATRIX%NUMBER_OF_COLUMNS
                    sum=sum+equations_matrix%ELEMENT_MATRIX%MATRIX(row_idx,column_idx)
                    equations_matrix%ELEMENT_MATRIX%MATRIX(row_idx,column_idx)=0.0_dp
                  ENDDO !column_idx
                  equations_matrix%ELEMENT_MATRIX%MATRIX(row_idx,row_idx)=sum
                  !Add the element matrice into the distributed equations matrix
                  CALL distributed_matrix_values_add(equations_matrix%MATRIX,equations_matrix%ELEMENT_MATRIX%ROW_DOFS(row_idx), &
                    & equations_matrix%ELEMENT_MATRIX%COLUMN_DOFS(row_idx),equations_matrix%ELEMENT_MATRIX%MATRIX(row_idx, &
                    & row_idx),err,error,*999)
                ENDDO !row_idx
              ELSE
                !Add the element matrice into the distributed equations matrix
                CALL distributed_matrix_values_add(equations_matrix%MATRIX,equations_matrix%ELEMENT_MATRIX%ROW_DOFS(1: &
                  & equations_matrix%ELEMENT_MATRIX%NUMBER_OF_ROWS),equations_matrix%ELEMENT_MATRIX%COLUMN_DOFS(1: &
                  & equations_matrix%ELEMENT_MATRIX%NUMBER_OF_COLUMNS),equations_matrix%ELEMENT_MATRIX%MATRIX(1: &
                  & equations_matrix%ELEMENT_MATRIX%NUMBER_OF_ROWS,1:equations_matrix%ELEMENT_MATRIX%NUMBER_OF_COLUMNS), &
                  & err,error,*999)
              ENDIF
            ENDIF
          ELSE
            local_error="Equations matrix for dynamic matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))// &
              & " is not associated."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ENDDO !matrix_idx
      ENDIF
      linear_matrices=>equations_matrices%LINEAR_MATRICES
      IF(ASSOCIATED(linear_matrices)) THEN
        !Add the element matrices
        DO matrix_idx=1,linear_matrices%NUMBER_OF_LINEAR_MATRICES
          equations_matrix=>linear_matrices%MATRICES(matrix_idx)%PTR
          IF(ASSOCIATED(equations_matrix)) THEN
            IF(equations_matrix%UPDATE_MATRIX) THEN
              !Handle lumped matrices
              IF(equations_matrix%LUMPED) THEN
                DO row_idx=1,equations_matrix%ELEMENT_MATRIX%NUMBER_OF_ROWS
                  sum=0.0_dp
                  DO column_idx=1,equations_matrix%ELEMENT_MATRIX%NUMBER_OF_COLUMNS
                    sum=sum+equations_matrix%ELEMENT_MATRIX%MATRIX(row_idx,column_idx)
                    equations_matrix%ELEMENT_MATRIX%MATRIX(row_idx,column_idx)=0.0_dp
                  ENDDO !column_idx
                  equations_matrix%ELEMENT_MATRIX%MATRIX(row_idx,row_idx)=sum
                  !Add the element matrice into the distributed equations matrix
                  CALL distributed_matrix_values_add(equations_matrix%MATRIX,equations_matrix%ELEMENT_MATRIX%ROW_DOFS(row_idx), &
                    & equations_matrix%ELEMENT_MATRIX%COLUMN_DOFS(row_idx),equations_matrix%ELEMENT_MATRIX%MATRIX(row_idx, &
                    & row_idx),err,error,*999)
                ENDDO !row_idx
              ELSE
                !Add the element matrice into the distributed equations matrix
                CALL distributed_matrix_values_add(equations_matrix%MATRIX,equations_matrix%ELEMENT_MATRIX%ROW_DOFS(1: &
                  & equations_matrix%ELEMENT_MATRIX%NUMBER_OF_ROWS),equations_matrix%ELEMENT_MATRIX%COLUMN_DOFS(1: &
                  & equations_matrix%ELEMENT_MATRIX%NUMBER_OF_COLUMNS),equations_matrix%ELEMENT_MATRIX%MATRIX(1: &
                  & equations_matrix%ELEMENT_MATRIX%NUMBER_OF_ROWS,1:equations_matrix%ELEMENT_MATRIX%NUMBER_OF_COLUMNS), &
                  & err,error,*999)
              ENDIF
            ENDIF
          ELSE
            local_error="Equations matrix for linear matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))// &
              & " is not associated."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ENDDO !matrix_idx
      ENDIF
      nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
      IF(ASSOCIATED(nonlinear_matrices)) THEN
        IF(nonlinear_matrices%UPDATE_RESIDUAL) THEN
          !Add the residual element vector
          CALL distributed_vector_values_add(nonlinear_matrices%RESIDUAL,nonlinear_matrices%ELEMENT_RESIDUAL%ROW_DOFS(1: &
            & nonlinear_matrices%ELEMENT_RESIDUAL%NUMBER_OF_ROWS),nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(1:nonlinear_matrices% &
            & element_residual%NUMBER_OF_ROWS),err,error,*999)
        ENDIF
      ENDIF
      rhs_vector=>equations_matrices%RHS_VECTOR
      IF(ASSOCIATED(rhs_vector)) THEN
        IF(rhs_vector%UPDATE_VECTOR) THEN
          !Add the rhs element vector
          CALL distributed_vector_values_add(rhs_vector%VECTOR,rhs_vector%ELEMENT_VECTOR%ROW_DOFS(1: &
            & rhs_vector%ELEMENT_VECTOR%NUMBER_OF_ROWS),rhs_vector%ELEMENT_VECTOR%VECTOR(1:rhs_vector% &
            & element_vector%NUMBER_OF_ROWS),err,error,*999)
        ENDIF
      ENDIF
      source_vector=>equations_matrices%SOURCE_VECTOR
      IF(ASSOCIATED(source_vector)) THEN
        IF(source_vector%UPDATE_VECTOR) THEN
          !Add the rhs element vector
          CALL distributed_vector_values_add(source_vector%VECTOR,source_vector%ELEMENT_VECTOR%ROW_DOFS(1: &
            & source_vector%ELEMENT_VECTOR%NUMBER_OF_ROWS),source_vector%ELEMENT_VECTOR%VECTOR(1:source_vector% &
            & element_vector%NUMBER_OF_ROWS),err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not allocated.",err,error,*999)
    ENDIF
#ifdef TAUPROF
    CALL tau_static_phase_stop("EQUATIONS_MATRICES_ELEMENT_ADD()")
#endif
    
    exits("EQUATIONS_MATRICES_ELEMENT_ADD")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_ELEMENT_ADD",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_element_add

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

  SUBROUTINE equations_matrices_element_calculate(EQUATIONS_MATRICES,ELEMENT_NUMBER,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
    TYPE(equations_jacobian_type), POINTER :: JACOBIAN_MATRIX
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_mapping_dynamic_type), POINTER :: DYNAMIC_MAPPING
    TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
    TYPE(equations_mapping_nonlinear_type), POINTER :: NONLINEAR_MAPPING
    TYPE(equations_mapping_rhs_type), POINTER :: RHS_MAPPING
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
    TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
    TYPE(equations_matrices_source_type), POINTER :: SOURCE_VECTOR
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE,COL_FIELD_VARIABLE
    TYPE(varying_string) :: LOCAL_ERROR

#ifdef TAUPROF
    CALL tau_static_phase_start("EQUATIONS_MATRICES_ELEMENT_CALCULATE()")
#endif

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

    IF(ASSOCIATED(equations_matrices)) THEN
      equations_mapping=>equations_matrices%EQUATIONS_MAPPING
      IF(ASSOCIATED(equations_mapping)) THEN
        dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
        IF(ASSOCIATED(dynamic_matrices)) THEN
          !Calculate the row and columns for the dynamic equations matrices
          dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
          IF(ASSOCIATED(dynamic_mapping)) THEN
            DO matrix_idx=1,dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES
              equations_matrix=>dynamic_matrices%MATRICES(matrix_idx)%PTR
              IF(ASSOCIATED(equations_matrix)) THEN
                field_variable=>dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrix_idx)%VARIABLE
                CALL equations_matrices_element_matrix_calculate(equations_matrix%ELEMENT_MATRIX,equations_matrix%UPDATE_MATRIX, &
                  & [element_number],[element_number],field_variable,field_variable,err,error,*999)
              ELSE
                local_error="Equations matrix for dynamic matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))// &
                  & " is not associated."
                CALL flagerror(local_error,err,error,*999)
              ENDIF
            ENDDO !matrix_idx
          ELSE
            CALL flagerror("Equations mapping dynamic mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        linear_matrices=>equations_matrices%LINEAR_MATRICES
        IF(ASSOCIATED(linear_matrices)) THEN
          !Calculate the row and columns for the linear equations matrices
          linear_mapping=>equations_mapping%LINEAR_MAPPING
          IF(ASSOCIATED(linear_mapping)) THEN
            DO matrix_idx=1,linear_matrices%NUMBER_OF_LINEAR_MATRICES
              equations_matrix=>linear_matrices%MATRICES(matrix_idx)%PTR
              IF(ASSOCIATED(equations_matrix)) THEN
                field_variable=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrix_idx)%VARIABLE
                CALL equations_matrices_element_matrix_calculate(equations_matrix%ELEMENT_MATRIX,equations_matrix%UPDATE_MATRIX, &
                  & [element_number],[element_number],field_variable,field_variable,err,error,*999)
              ELSE
                local_error="Equations matrix for linear matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))// &
                  & " is not associated."
                CALL flagerror(local_error,err,error,*999)
              ENDIF
            ENDDO !matrix_idx
          ELSE
            CALL flagerror("Equations mapping linear mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
        IF(ASSOCIATED(nonlinear_matrices)) THEN
          !Calculate the rows and columns of the Jacobian
          nonlinear_mapping=>equations_mapping%NONLINEAR_MAPPING
          IF(ASSOCIATED(nonlinear_mapping)) THEN
            field_variable=>nonlinear_mapping%JACOBIAN_TO_VAR_MAP(1)%VARIABLE !Row field variable
            DO matrix_idx=1,nonlinear_matrices%NUMBER_OF_JACOBIANS
              jacobian_matrix=>nonlinear_matrices%JACOBIANS(matrix_idx)%PTR
              IF(ASSOCIATED(jacobian_matrix)) THEN
                col_field_variable=>nonlinear_mapping%JACOBIAN_TO_VAR_MAP(matrix_idx)%VARIABLE
                CALL equations_matrices_element_matrix_calculate(jacobian_matrix%ELEMENT_JACOBIAN,jacobian_matrix%UPDATE_JACOBIAN, &
                  & [element_number],[element_number],field_variable,col_field_variable,err,error,*999)
              ELSE
                CALL flagerror("Jacobian matrix is not associated.",err,error,*999)
              ENDIF
            ENDDO
            !Calculate the rows of the equations residual
            rhs_mapping=>equations_mapping%RHS_MAPPING
            IF(ASSOCIATED(rhs_mapping)) THEN
              field_variable=>rhs_mapping%RHS_VARIABLE
            ELSE
              field_variable=>nonlinear_mapping%JACOBIAN_TO_VAR_MAP(1)%VARIABLE
            ENDIF
            CALL equations_matrices_element_vector_calculate(nonlinear_matrices%ELEMENT_RESIDUAL,nonlinear_matrices% &
              & update_residual,element_number,field_variable,err,error,*999)
            nonlinear_matrices%ELEMENT_RESIDUAL_CALCULATED=0
          ELSE
            CALL flagerror("Equations mapping nonlinear mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        rhs_vector=>equations_matrices%RHS_VECTOR
        IF(ASSOCIATED(rhs_vector)) THEN
          rhs_mapping=>equations_mapping%RHS_MAPPING
          IF(ASSOCIATED(rhs_mapping)) THEN
            !Calculate the rows  for the equations RHS
            field_variable=>rhs_mapping%RHS_VARIABLE
            CALL equations_matrices_element_vector_calculate(rhs_vector%ELEMENT_VECTOR,rhs_vector%UPDATE_VECTOR,element_number, &
              & field_variable,err,error,*999)
          ELSE
            CALL flagerror("Equations mapping rhs mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        source_vector=>equations_matrices%SOURCE_VECTOR
        IF(ASSOCIATED(source_vector)) THEN
          !Calculate the rows the equations source. The number of rows is not set by the source field so take the number of rows
          !from the RHS vector in the first instance.
          rhs_mapping=>equations_mapping%RHS_MAPPING
          IF(ASSOCIATED(rhs_mapping)) THEN
            field_variable=>rhs_mapping%RHS_VARIABLE
            CALL equations_matrices_element_vector_calculate(source_vector%ELEMENT_VECTOR,source_vector%UPDATE_VECTOR, &
              & element_number,field_variable,err,error,*999)
          ELSE
            CALL flagerror("Equations mapping rhs mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
      ELSE
        CALL flagerror("Equations mapping is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not allocated",err,error,*999)
    ENDIF
    
#ifdef TAUPROF
    CALL tau_static_phase_stop("EQUATIONS_MATRICES_ELEMENT_CALCULATE()")
#endif
    
    exits("EQUATIONS_MATRICES_ELEMENT_CALCULATE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_ELEMENT_CALCULATE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_element_calculate

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

  SUBROUTINE equationsmatrices_nodalcalculate(equationsMatrices,nodeNumber,err,error,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: equationsMatrices
    INTEGER(INTG), INTENT(IN) :: nodeNumber
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: matrixIdx
    TYPE(equations_jacobian_type), POINTER :: jacobianMatrix
    TYPE(equations_mapping_type), POINTER :: equationsMapping
    TYPE(equations_mapping_dynamic_type), POINTER :: dynamicMapping
    TYPE(equations_mapping_linear_type), POINTER :: linearMapping
    TYPE(equations_mapping_nonlinear_type), POINTER :: nonlinearMapping
    TYPE(equations_mapping_rhs_type), POINTER :: rhsMapping
    TYPE(equations_matrices_dynamic_type), POINTER :: dynamicMatrices
    TYPE(equations_matrices_linear_type), POINTER :: linearMatrices
    TYPE(equations_matrices_nonlinear_type), POINTER :: nonlinearMatrices
    TYPE(equations_matrices_rhs_type), POINTER :: rhsVector
    TYPE(equations_matrices_source_type), POINTER :: sourceVector
    TYPE(equations_matrix_type), POINTER :: equationsMatrix
    TYPE(field_variable_type), POINTER :: fieldVariable,columnFieldVariable
    TYPE(varying_string) :: localError

#ifdef TAUPROF
    CALL tau_static_phase_start("EquationsMatrices_NodalCalculate()")
#endif

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

    IF(ASSOCIATED(equationsmatrices)) THEN
      equationsmapping=>equationsmatrices%EQUATIONS_MAPPING
      IF(ASSOCIATED(equationsmapping)) THEN
        dynamicmatrices=>equationsmatrices%DYNAMIC_MATRICES
        IF(ASSOCIATED(dynamicmatrices)) THEN
          !Calculate the row and columns for the dynamic equations matrices
          dynamicmapping=>equationsmapping%DYNAMIC_MAPPING
          IF(ASSOCIATED(dynamicmapping)) THEN
            DO matrixidx=1,dynamicmatrices%NUMBER_OF_DYNAMIC_MATRICES
              equationsmatrix=>dynamicmatrices%MATRICES(matrixidx)%PTR
              IF(ASSOCIATED(equationsmatrix)) THEN
                fieldvariable=>dynamicmapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrixidx)%VARIABLE
                CALL equationsmatrices_nodalmatrixcalculate(equationsmatrix%NodalMatrix,equationsmatrix%UPDATE_MATRIX, &
                  & nodenumber,nodenumber,fieldvariable,fieldvariable,err,error,*999)
              ELSE
                localerror="Equations matrix for dynamic matrix number "//trim(numbertovstring(matrixidx,"*",err,error))// &
                  & " is not associated."
                CALL flagerror(localerror,err,error,*999)
              ENDIF
            ENDDO !matrixIdx
          ELSE
            CALL flagerror("Equations mapping dynamic mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        linearmatrices=>equationsmatrices%LINEAR_MATRICES
        IF(ASSOCIATED(linearmatrices)) THEN
          !Calculate the row and columns for the linear equations matrices
          linearmapping=>equationsmapping%LINEAR_MAPPING
          IF(ASSOCIATED(linearmapping)) THEN
            DO matrixidx=1,linearmatrices%NUMBER_OF_LINEAR_MATRICES
              equationsmatrix=>linearmatrices%MATRICES(matrixidx)%PTR
              IF(ASSOCIATED(equationsmatrix)) THEN
                fieldvariable=>linearmapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrixidx)%VARIABLE
                CALL equationsmatrices_nodalmatrixcalculate(equationsmatrix%NodalMatrix,equationsmatrix%UPDATE_MATRIX, &
                  & nodenumber,nodenumber,fieldvariable,fieldvariable,err,error,*999)
              ELSE
                localerror="Equations matrix for linear matrix number "//trim(numbertovstring(matrixidx,"*",err,error))// &
                  & " is not associated."
                CALL flagerror(localerror,err,error,*999)
              ENDIF
            ENDDO !matrixIdx
          ELSE
            CALL flagerror("Equations mapping linear mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        nonlinearmatrices=>equationsmatrices%NONLINEAR_MATRICES
        IF(ASSOCIATED(nonlinearmatrices)) THEN
          !Calculate the rows and columns of the Jacobian
          nonlinearmapping=>equationsmapping%NONLINEAR_MAPPING
          IF(ASSOCIATED(nonlinearmapping)) THEN
            fieldvariable=>nonlinearmapping%JACOBIAN_TO_VAR_MAP(1)%VARIABLE !Row field variable
            DO matrixidx=1,nonlinearmatrices%NUMBER_OF_JACOBIANS
              jacobianmatrix=>nonlinearmatrices%JACOBIANS(matrixidx)%PTR
              IF(ASSOCIATED(jacobianmatrix)) THEN
                columnfieldvariable=>nonlinearmapping%JACOBIAN_TO_VAR_MAP(matrixidx)%VARIABLE
                CALL equationsmatrices_nodalmatrixcalculate(jacobianmatrix%NodalJacobian,jacobianmatrix%UPDATE_JACOBIAN, &
                  & nodenumber,nodenumber,fieldvariable,columnfieldvariable,err,error,*999)
              ELSE
                CALL flagerror("Jacobian matrix is not associated.",err,error,*999)
              ENDIF
            ENDDO
            !Calculate the rows of the equations residual
            rhsmapping=>equationsmapping%RHS_MAPPING
            IF(ASSOCIATED(rhsmapping)) THEN
              fieldvariable=>rhsmapping%RHS_VARIABLE
            ELSE
              fieldvariable=>nonlinearmapping%JACOBIAN_TO_VAR_MAP(1)%VARIABLE
            ENDIF
            CALL equationsmatrices_nodalvectorcalculate(nonlinearmatrices%NodalResidual,nonlinearmatrices% &
              & update_residual,nodenumber,fieldvariable,err,error,*999)
            nonlinearmatrices%NodalResidualCalculated=0
          ELSE
            CALL flagerror("Equations mapping nonlinear mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        rhsvector=>equationsmatrices%RHS_VECTOR
        IF(ASSOCIATED(rhsvector)) THEN
          rhsmapping=>equationsmapping%RHS_MAPPING
          IF(ASSOCIATED(rhsmapping)) THEN
            !Calculate the rows  for the equations RHS
            fieldvariable=>rhsmapping%RHS_VARIABLE
            CALL equationsmatrices_nodalvectorcalculate(rhsvector%NodalVector,rhsvector%UPDATE_VECTOR,nodenumber, &
              & fieldvariable,err,error,*999)
          ELSE
            CALL flagerror("Equations mapping rhs mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        sourcevector=>equationsmatrices%SOURCE_VECTOR
        IF(ASSOCIATED(sourcevector)) THEN
          !Calculate the rows the equations source. The number of rows is not set by the source field so take the number of rows
          !from the RHS vector in the first instance.
          rhsmapping=>equationsmapping%RHS_MAPPING
          IF(ASSOCIATED(rhsmapping)) THEN
            fieldvariable=>rhsmapping%RHS_VARIABLE
            CALL equationsmatrices_nodalvectorcalculate(sourcevector%NodalVector,sourcevector%UPDATE_VECTOR, &
              & nodenumber,fieldvariable,err,error,*999)
          ELSE
            CALL flagerror("Equations mapping rhs mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
      ELSE
        CALL flagerror("Equations mapping is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not allocated",err,error,*999)
    ENDIF
    
#ifdef TAUPROF
    CALL tau_static_phase_stop("EquationsMatrices_NodalCalculate()")
#endif
    
    exits("EquationsMatrices_NodalCalculate")
    RETURN
999 errorsexits("EquationsMatrices_NodalCalculate",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_nodalcalculate

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

  SUBROUTINE equationsmatrices_nodalmatrixcalculate(nodalMatrix,updateMatrix,rowNodeNumber,columnNodeNumber, &
    & rowsfieldvariable,colsfieldvariable,err,error,*)

    !Argument variables
    TYPE(nodalmatrixtype) :: nodalMatrix
    LOGICAL :: updateMatrix
    INTEGER(INTG), INTENT(IN) :: rowNodeNumber
    INTEGER(INTG), INTENT(IN) :: columnNodeNumber
    TYPE(field_variable_type), POINTER :: rowsFieldVariable
    TYPE(field_variable_type), POINTER :: colsFieldVariable
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: componentIdx
    INTEGER(INTG) :: localRow,globalRow,localColumn,globalColumn
    INTEGER(INTG) :: numberOfDerivatives,numberOfVersions,versionIdx,derivativeIdx
    TYPE(domain_nodes_type), POINTER :: nodesTopology
    TYPE(varying_string) :: localError

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

    IF(ASSOCIATED(rowsfieldvariable)) THEN
      IF(ASSOCIATED(colsfieldvariable)) THEN
        nodalmatrix%numberOfRows=0
        nodalmatrix%numberOfColumns=0
        IF(updatematrix) THEN
          IF(ASSOCIATED(rowsfieldvariable,colsfieldvariable)) THEN
            !Row and columns variable is the same.
            DO componentidx=1,rowsfieldvariable%NUMBER_OF_COMPONENTS
              nodestopology=>rowsfieldvariable%COMPONENTS(componentidx)%DOMAIN%TOPOLOGY%NODES
              IF(rownodenumber>=1.AND.rownodenumber<=nodestopology%TOTAL_NUMBER_OF_NODES) THEN
                SELECT CASE(rowsfieldvariable%COMPONENTS(componentidx)%INTERPOLATION_TYPE)
                CASE(field_constant_interpolation)
                  localrow=rowsfieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
                  globalrow=rowsfieldvariable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(localrow)
                  nodalmatrix%numberOfRows=nodalmatrix%numberOfRows+1
                  nodalmatrix%numberOfColumns=nodalmatrix%numberOfColumns+1
                  nodalmatrix%rowDofs(nodalmatrix%numberOfRows)=localrow
                  nodalmatrix%columnDofs(nodalmatrix%numberOfColumns)=globalrow
                CASE(field_element_based_interpolation)
                  localrow=rowsfieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP% &
                    & elements(rownodenumber)
                  globalrow=rowsfieldvariable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(localrow)
                  nodalmatrix%numberOfRows=nodalmatrix%numberOfRows+1
                  nodalmatrix%numberOfColumns=nodalmatrix%numberOfColumns+1
                  nodalmatrix%rowDofs(nodalmatrix%numberOfRows)=localrow
                  nodalmatrix%columnDofs(nodalmatrix%numberOfColumns)=globalrow
                CASE(field_node_based_interpolation)
                  numberofderivatives=rowsfieldvariable%components(componentidx)%domain%topology%nodes%nodes(rownodenumber)% &
                    & number_of_derivatives
                  DO derivativeidx=1,numberofderivatives
                    numberofversions=rowsfieldvariable%COMPONENTS(componentidx)%DOMAIN%TOPOLOGY%NODES%NODES(rownodenumber)% &
                      & derivatives(derivativeidx)%numberOfVersions
                    DO versionidx=1,numberofversions
                      localrow=rowsfieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% & 
                        & nodes(rownodenumber)%DERIVATIVES(derivativeidx)%VERSIONS(versionidx)
                      globalrow=rowsfieldvariable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(localrow)
                      nodalmatrix%numberOfRows=nodalmatrix%numberOfRows+1
                      nodalmatrix%numberOfColumns=nodalmatrix%numberOfColumns+1
                      nodalmatrix%rowDofs(nodalmatrix%numberOfRows)=localrow
                      nodalmatrix%columnDofs(nodalmatrix%numberOfColumns)=globalrow
                    ENDDO !versionIdx
                  ENDDO !derivativeIdx
                CASE(field_grid_point_based_interpolation)
                  CALL flagerror("Not implemented.",err,error,*999)
                CASE(field_gauss_point_based_interpolation)
                  CALL flagerror("Not implemented.",err,error,*999)
                CASE DEFAULT
                  localerror="The interpolation type of "// &
                    & trim(numbertovstring(rowsfieldvariable%COMPONENTS(componentidx)%INTERPOLATION_TYPE,"*",err,error))// &
                    & " is invalid for component number "// &
                    & trim(numbertovstring(componentidx,"*",err,error))// &
                    & " of rows field variable type "// &
                    & trim(numbertovstring(rowsfieldvariable%VARIABLE_TYPE,"*",err,error))//"."
                  CALL flagerror(localerror,err,error,*999)          
                END SELECT
              ELSE
                localerror="Nodal number "//trim(numbertovstring(rownodenumber,"*",err,error))// &
                  & " is invalid for component number "//trim(numbertovstring(componentidx,"*",err,error))// &
                  & " of rows field variable type "// &
                  & trim(numbertovstring(rowsfieldvariable%VARIABLE_TYPE,"*",err,error))// &
                  & ". The nodal number must be between 1 and "// &
                  & trim(numbertovstring(nodestopology%TOTAL_NUMBER_OF_NODES,"*",err,error))//"."
                CALL flagerror(localerror,err,error,*999)
              ENDIF
            ENDDO !componentIdx
          ELSE
            !Row and column variables are different
            !Row mapping
            DO componentidx=1,rowsfieldvariable%NUMBER_OF_COMPONENTS
              nodestopology=>rowsfieldvariable%COMPONENTS(componentidx)%DOMAIN%TOPOLOGY%NODES
              IF(rownodenumber>=1.AND.rownodenumber<=nodestopology%TOTAL_NUMBER_OF_NODES) THEN
                SELECT CASE(rowsfieldvariable%COMPONENTS(componentidx)%INTERPOLATION_TYPE)
                CASE(field_constant_interpolation)
                  localrow=rowsfieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
                  nodalmatrix%numberOfRows=nodalmatrix%numberOfRows+1
                  nodalmatrix%rowDofs(nodalmatrix%numberOfRows)=localrow
                CASE(field_element_based_interpolation)
                  localrow=rowsfieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP% &
                    & elements(rownodenumber)
                  nodalmatrix%numberOfRows=nodalmatrix%numberOfRows+1
                  nodalmatrix%rowDofs(nodalmatrix%numberOfRows)=localrow
                CASE(field_node_based_interpolation)
                  numberofderivatives=rowsfieldvariable%components(componentidx)%domain%topology%nodes%nodes(rownodenumber)% &
                    & number_of_derivatives
                  DO derivativeidx=1,numberofderivatives
                    numberofversions=colsfieldvariable%COMPONENTS(componentidx)%DOMAIN%TOPOLOGY%NODES%NODES(rownodenumber)% &
                      & derivatives(derivativeidx)%numberOfVersions
                    DO versionidx=1,numberofversions
                      localrow=rowsfieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% & 
                        & nodes(rownodenumber)%DERIVATIVES(derivativeidx)%VERSIONS(versionidx)
                      nodalmatrix%numberOfRows=nodalmatrix%numberOfRows+1
                      nodalmatrix%rowDofs(nodalmatrix%numberOfRows)=localrow
                    ENDDO !versionIdx
                  ENDDO !derivativeIdx
                CASE(field_grid_point_based_interpolation)
                  CALL flagerror("Not implemented.",err,error,*999)
                CASE(field_gauss_point_based_interpolation)
                  CALL flagerror("Not implemented.",err,error,*999)
                CASE DEFAULT
                  localerror="The interpolation type of "// &
                    & trim(numbertovstring(rowsfieldvariable%COMPONENTS(componentidx)%INTERPOLATION_TYPE,"*",err,error))// &
                    & " is invalid for component number "// &
                    & trim(numbertovstring(componentidx,"*",err,error))// &
                    & " of rows field variable type "// &
                    & trim(numbertovstring(rowsfieldvariable%VARIABLE_TYPE,"*",err,error))//"."
                  CALL flagerror(localerror,err,error,*999)          
                END SELECT
              ELSE
                localerror="Row nodal number "//trim(numbertovstring(rownodenumber,"*",err,error))// &
                  & " is invalid for component number "//trim(numbertovstring(componentidx,"*",err,error))// &
                  & " of rows field variable type "// &
                  & trim(numbertovstring(rowsfieldvariable%VARIABLE_TYPE,"*",err,error))// &
                  & ". The nodal number must be between 1 and "// &
                  & trim(numbertovstring(nodestopology%TOTAL_NUMBER_OF_NODES,"*",err,error))//"."
                CALL flagerror(localerror,err,error,*999)
              ENDIF
            ENDDO !componentIdx
            !Column mapping
            DO componentidx=1,colsfieldvariable%NUMBER_OF_COMPONENTS
              nodestopology=>colsfieldvariable%COMPONENTS(componentidx)%DOMAIN%TOPOLOGY%NODES
              IF(columnnodenumber>=1.AND.columnnodenumber<=nodestopology%TOTAL_NUMBER_OF_NODES) THEN
                SELECT CASE(colsfieldvariable%COMPONENTS(componentidx)%INTERPOLATION_TYPE)
                CASE(field_constant_interpolation)
                  localcolumn=colsfieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
                  globalcolumn=colsfieldvariable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(localcolumn)
                  nodalmatrix%numberOfColumns=nodalmatrix%numberOfColumns+1
                  nodalmatrix%columnDofs(nodalmatrix%numberOfColumns)=globalcolumn
                CASE(field_element_based_interpolation)
                  localcolumn=colsfieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP% &
                    & elements(columnnodenumber)
                  globalcolumn=colsfieldvariable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(localcolumn)
                  nodalmatrix%numberOfColumns=nodalmatrix%numberOfColumns+1
                  nodalmatrix%columnDofs(nodalmatrix%numberOfColumns)=globalcolumn
                CASE(field_node_based_interpolation)
                  numberofderivatives=colsfieldvariable%components(componentidx)%domain%topology%nodes%nodes(rownodenumber)% &
                    & number_of_derivatives
                  DO derivativeidx=1,numberofderivatives
                    numberofversions=colsfieldvariable%COMPONENTS(componentidx)%DOMAIN%TOPOLOGY%NODES%NODES(rownodenumber)% &
                      & derivatives(derivativeidx)%numberOfVersions
                    DO versionidx=1,numberofversions
                      localrow=colsfieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% & 
                        & nodes(rownodenumber)%DERIVATIVES(derivativeidx)%VERSIONS(versionidx)
                      nodalmatrix%numberOfColumns=nodalmatrix%numberOfColumns+1
                      nodalmatrix%columnDofs(nodalmatrix%numberOfColumns)=localrow
                    ENDDO !versionIdx
                  ENDDO !derivativeIdx
                CASE(field_grid_point_based_interpolation)
                  CALL flagerror("Not implemented.",err,error,*999)
                CASE(field_gauss_point_based_interpolation)
                  CALL flagerror("Not implemented.",err,error,*999)
                CASE DEFAULT
                  localerror="The interpolation type of "// &
                    & trim(numbertovstring(colsfieldvariable%COMPONENTS(componentidx)%INTERPOLATION_TYPE,"*",err,error))// &
                    & " is invalid for component number "// &
                    & trim(numbertovstring(componentidx,"*",err,error))// &
                    & " of column field variable type "// &
                    & trim(numbertovstring(colsfieldvariable%VARIABLE_TYPE,"*",err,error))//"."
                  CALL flagerror(localerror,err,error,*999)          
                END SELECT
              ELSE
                localerror="Column nodal number "//trim(numbertovstring(columnnodenumber,"*",err,error))// &
                  & " is invalid for component number "//trim(numbertovstring(componentidx,"*",err,error))// &
                  & " of column field variable type "// &
                  & trim(numbertovstring(colsfieldvariable%VARIABLE_TYPE,"*",err,error))// &
                  & ". The nodal number must be between 1 and "// &
                  & trim(numbertovstring(nodestopology%TOTAL_NUMBER_OF_NODES,"*",err,error))//"."
                CALL flagerror(localerror,err,error,*999)
              ENDIF
            ENDDO !componentIdx
          ENDIF
          nodalmatrix%matrix=0.0_dp
        ENDIF
      ELSE
        CALL flagerror("Columns field variable is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Rows field variable is not associated.",err,error,*999)
    ENDIF
    
    exits("EquationsMatrices_NodalMatrixCalculate")
    RETURN
999 errorsexits("EquationsMatrices_NodalMatrixCalculate",err,error)
    RETURN 1
    
  END SUBROUTINE equationsmatrices_nodalmatrixcalculate

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

  SUBROUTINE equationsmatrices_nodalvectorcalculate(nodalVector,updateVector,rowNodeNumber,rowsFieldVariable, &
    & err,error,*)

    !Argument variables
    TYPE(nodalvectortype) :: nodalVector
    LOGICAL :: updateVector
    INTEGER(INTG), INTENT(IN) :: rowNodeNumber
    TYPE(field_variable_type), POINTER :: rowsFieldVariable
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: componentIdx,localRow
    INTEGER(INTG) :: numberOfDerivatives,numberOfVersions,versionIdx,derivativeIdx
    TYPE(domain_nodes_type), POINTER :: nodesTopology
    TYPE(varying_string) :: localError
    
    enters("EquationsMatrices_NodalVectorCalculate",err,error,*999)

    IF(ASSOCIATED(rowsfieldvariable)) THEN
      !Calculate the rows for the nodal vector
      nodalvector%numberOfRows=0
      IF(updatevector) THEN
        DO componentidx=1,rowsfieldvariable%NUMBER_OF_COMPONENTS
          nodestopology=>rowsfieldvariable%components(componentidx)%domain%topology%nodes
          IF(rownodenumber>=1.AND.rownodenumber<=nodestopology%TOTAL_NUMBER_OF_NODES) THEN
            SELECT CASE(rowsfieldvariable%COMPONENTS(componentidx)%INTERPOLATION_TYPE)
            CASE(field_constant_interpolation)
              localrow=rowsfieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
              nodalvector%numberOfRows=nodalvector%numberOfRows+1
              nodalvector%rowDofs(nodalvector%numberOfRows)=localrow
            CASE(field_element_based_interpolation)
              localrow=rowsfieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP% &
                & elements(rownodenumber)
              nodalvector%numberOfRows=nodalvector%numberOfRows+1
              nodalvector%rowDofs(nodalvector%numberOfRows)=localrow
            CASE(field_node_based_interpolation)
              numberofderivatives=rowsfieldvariable%components(componentidx)%domain%topology%nodes%nodes(rownodenumber)% &
                & number_of_derivatives
              DO derivativeidx=1,numberofderivatives
                numberofversions=rowsfieldvariable%COMPONENTS(componentidx)%DOMAIN%TOPOLOGY%NODES%NODES(rownodenumber)% &
                  & derivatives(derivativeidx)%numberOfVersions
                DO versionidx=1,numberofversions
                  localrow=rowsfieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% & 
                    & nodes(rownodenumber)%derivatives(derivativeidx)%versions(versionidx)
                  nodalvector%numberOfRows=nodalvector%numberOfRows+1
                  nodalvector%rowDofs(nodalvector%numberOfRows)=localrow
                ENDDO !versionIdx
              ENDDO !derivativeIdx
            CASE(field_grid_point_based_interpolation)
              CALL flagerror("Not implemented.",err,error,*999)
            CASE(field_gauss_point_based_interpolation)
              CALL flagerror("Not implemented.",err,error,*999)
            CASE DEFAULT
              localerror="The interpolation type of "// &
                & trim(numbertovstring(rowsfieldvariable%COMPONENTS(componentidx)%INTERPOLATION_TYPE,"*",err,error))// &
                & " is invalid for component number "// &
                & trim(numbertovstring(componentidx,"*",err,error))// &
                & " of rows field variable type "// &
                & trim(numbertovstring(rowsfieldvariable%VARIABLE_TYPE,"*",err,error))//"."
              CALL flagerror(localerror,err,error,*999)          
            END SELECT
          ELSE
            localerror="Node number "//trim(numbertovstring(rownodenumber,"*",err,error))// &
              & " is invalid for component number "//trim(numbertovstring(componentidx,"*",err,error))// &
              & " of rows field variable type "//trim(numbertovstring(rowsfieldvariable%VARIABLE_TYPE,"*",err,error))// &
              & ". The nodal number must be between 1 and "// &
              & trim(numbertovstring(nodestopology%TOTAL_NUMBER_OF_NODES,"*",err,error))//"."
            CALL flagerror(localerror,err,error,*999)
          ENDIF
        ENDDO !componentIdx
        nodalvector%vector=0.0_dp
      ENDIF
    ELSE
      CALL flagerror("Rows field variable is not associated.",err,error,*999)
    ENDIF
    
    exits("EquationsMatrices_NodalVectorCalculate")
    RETURN
999 errorsexits("EquationsMatrices_NodalVectorCalculate",err,error)
    RETURN 1
    
  END SUBROUTINE equationsmatrices_nodalvectorcalculate

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

  SUBROUTINE equationsmatrices_nodeadd(equationsMatrices,err,error,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: equationsMatrices
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: columnIdx,matrixIdx,rowIdx
    REAL(DP) :: sum
    TYPE(equations_matrices_dynamic_type), POINTER :: dynamicMatrices
    TYPE(equations_matrices_linear_type), POINTER :: linearMatrices
    TYPE(equations_matrices_nonlinear_type), POINTER :: nonlinearMatrices
    TYPE(equations_matrices_rhs_type), POINTER :: rhsVector
    TYPE(equations_matrices_source_type), POINTER :: sourceVector
    TYPE(equations_matrix_type), POINTER :: equationsMatrix
    TYPE(varying_string) :: localError
    
#ifdef TAUPROF
    CALL tau_static_phase_start("EquationsMatrices_NodeAdd()")
#endif

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

    IF(ASSOCIATED(equationsmatrices)) THEN
      dynamicmatrices=>equationsmatrices%DYNAMIC_MATRICES
      IF(ASSOCIATED(dynamicmatrices)) THEN
        !Add the nodal matrices
        DO matrixidx=1,dynamicmatrices%NUMBER_OF_DYNAMIC_MATRICES
          equationsmatrix=>dynamicmatrices%MATRICES(matrixidx)%PTR
          IF(ASSOCIATED(equationsmatrix)) THEN
            IF(equationsmatrix%UPDATE_MATRIX) THEN
              !Handle lumped matrices
              IF(equationsmatrix%LUMPED) THEN
                DO rowidx=1,equationsmatrix%NodalMatrix%numberOfRows
                  sum=0.0_dp
                  DO columnidx=1,equationsmatrix%NodalMatrix%numberOfColumns
                    sum=sum+equationsmatrix%NodalMatrix%matrix(rowidx,columnidx)
                    equationsmatrix%NodalMatrix%matrix(rowidx,columnidx)=0.0_dp
                  ENDDO !columnIdx
                  equationsmatrix%NodalMatrix%matrix(rowidx,rowidx)=sum
                  !Add the nodal matrice into the distributed equations matrix
                  CALL distributed_matrix_values_add(equationsmatrix%matrix,equationsmatrix%NodalMatrix%rowDofs(rowidx), &
                    & equationsmatrix%NodalMatrix%columnDofs(rowidx),equationsmatrix%NodalMatrix%matrix(rowidx, &
                    & rowidx),err,error,*999)
                ENDDO !rowIdx
              ELSE
                !Add the nodal matrice into the distributed equations matrix
                CALL distributed_matrix_values_add(equationsmatrix%matrix,equationsmatrix%NodalMatrix%rowDofs(1: &
                  & equationsmatrix%NodalMatrix%numberOfRows),equationsmatrix%NodalMatrix%columnDofs(1: &
                  & equationsmatrix%NodalMatrix%numberOfColumns),equationsmatrix%NodalMatrix%matrix(1: &
                  & equationsmatrix%NodalMatrix%numberOfRows,1:equationsmatrix%NodalMatrix%numberOfColumns), &
                  & err,error,*999)
              ENDIF
            ENDIF
          ELSE
            localerror="Equations matrix for dynamic matrix number "//trim(numbertovstring(matrixidx,"*",err,error))// &
              & " is not associated."
            CALL flagerror(localerror,err,error,*999)
          ENDIF
        ENDDO !matrixIdx
      ENDIF
      linearmatrices=>equationsmatrices%LINEAR_MATRICES
      IF(ASSOCIATED(linearmatrices)) THEN
        !Add the nodal matrices
        DO matrixidx=1,linearmatrices%NUMBER_OF_LINEAR_MATRICES
          equationsmatrix=>linearmatrices%MATRICES(matrixidx)%PTR
          IF(ASSOCIATED(equationsmatrix)) THEN
            IF(equationsmatrix%UPDATE_MATRIX) THEN
              !Handle lumped matrices
              IF(equationsmatrix%LUMPED) THEN
                DO rowidx=1,equationsmatrix%NodalMatrix%numberOfRows
                  sum=0.0_dp
                  DO columnidx=1,equationsmatrix%NodalMatrix%numberOfColumns
                    sum=sum+equationsmatrix%NodalMatrix%matrix(rowidx,columnidx)
                    equationsmatrix%NodalMatrix%matrix(rowidx,columnidx)=0.0_dp
                  ENDDO !columnIdx
                  equationsmatrix%NodalMatrix%matrix(rowidx,rowidx)=sum
                  !Add the nodal matrice into the distributed equations matrix
                  CALL distributed_matrix_values_add(equationsmatrix%matrix,equationsmatrix%NodalMatrix%rowDofs(rowidx), &
                    & equationsmatrix%NodalMatrix%columnDofs(rowidx),equationsmatrix%NodalMatrix%matrix(rowidx, &
                    & rowidx),err,error,*999)
                ENDDO !rowIdx
              ELSE
                !Add the nodal matrice into the distributed equations matrix
                CALL distributed_matrix_values_add(equationsmatrix%matrix,equationsmatrix%NodalMatrix%rowDofs(1: &
                  & equationsmatrix%NodalMatrix%numberOfRows),equationsmatrix%NodalMatrix%columnDofs(1: &
                  & equationsmatrix%NodalMatrix%numberOfColumns),equationsmatrix%NodalMatrix%matrix(1: &
                  & equationsmatrix%NodalMatrix%numberOfRows,1:equationsmatrix%NodalMatrix%numberOfColumns), &
                  & err,error,*999)
              ENDIF
            ENDIF
          ELSE
            localerror="Equations matrix for linear matrix number "//trim(numbertovstring(matrixidx,"*",err,error))// &
              & " is not associated."
            CALL flagerror(localerror,err,error,*999)
          ENDIF
        ENDDO !matrixIdx
      ENDIF
      nonlinearmatrices=>equationsmatrices%NONLINEAR_MATRICES
      IF(ASSOCIATED(nonlinearmatrices)) THEN
        IF(nonlinearmatrices%UPDATE_RESIDUAL) THEN
          !Add the residual nodal vector
          CALL distributed_vector_values_add(nonlinearmatrices%RESIDUAL,nonlinearmatrices%NodalResidual%rowDofs(1: &
            & nonlinearmatrices%NodalResidual%numberOfRows),nonlinearmatrices%NodalResidual%vector(1:nonlinearmatrices% &
            & nodalresidual%numberOfRows),err,error,*999)
        ENDIF
      ENDIF
      rhsvector=>equationsmatrices%RHS_VECTOR
      IF(ASSOCIATED(rhsvector)) THEN
        IF(rhsvector%UPDATE_VECTOR) THEN
          !Add the rhs nodal vector
          CALL distributed_vector_values_add(rhsvector%vector,rhsvector%NodalVector%rowDofs(1: &
            & rhsvector%NodalVector%numberOfRows),rhsvector%NodalVector%vector(1:rhsvector% &
            & nodalvector%numberOfRows),err,error,*999)
        ENDIF
      ENDIF
      sourcevector=>equationsmatrices%SOURCE_VECTOR
      IF(ASSOCIATED(sourcevector)) THEN
        IF(sourcevector%UPDATE_VECTOR) THEN
          !Add the rhs nodal vector
          CALL distributed_vector_values_add(sourcevector%vector,sourcevector%NodalVector%rowDofs(1: &
            & sourcevector%NodalVector%numberOfRows),sourcevector%NodalVector%vector(1:sourcevector% &
            & nodalvector%numberOfRows),err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not allocated.",err,error,*999)
    ENDIF
#ifdef TAUPROF
    CALL tau_static_phase_stop("EquationsMatrices_NodeAdd()")
#endif
    
    exits("EquationsMatrices_NodeAdd")
    RETURN
999 errorsexits("EquationsMatrices_NodeAdd",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_nodeadd

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

  SUBROUTINE equationsmatrices_nodalinitialise(equationsMatrices,err,error,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: equationsMatrices !The equations matrices to initialise the nodal information for
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: matrixIdx
    TYPE(equations_jacobian_type), POINTER :: jacobianMatrix
    TYPE(equations_mapping_type), POINTER :: equationsMapping
    TYPE(equations_mapping_dynamic_type), POINTER :: dynamicMapping
    TYPE(equations_mapping_linear_type), POINTER :: linearMapping
    TYPE(equations_mapping_nonlinear_type), POINTER :: nonlinearMapping
    TYPE(equations_mapping_rhs_type), POINTER :: rhsMapping
    TYPE(equations_matrices_dynamic_type), POINTER :: dynamicMatrices
    TYPE(equations_matrices_linear_type), POINTER :: linearMatrices
    TYPE(equations_matrices_nonlinear_type), POINTER :: nonlinearMatrices
    TYPE(equations_matrices_rhs_type), POINTER :: rhsVector
    TYPE(equations_matrices_source_type), POINTER :: sourceVector
    TYPE(equations_matrix_type), POINTER :: equationsMatrix
    TYPE(field_variable_type), POINTER :: fieldVariable,columnFieldVariable
    TYPE(varying_string) :: localError
    
    enters("EquationsMatrices_NodalInitialise",err,error,*999)

    IF(ASSOCIATED(equationsmatrices)) THEN
      equationsmapping=>equationsmatrices%EQUATIONS_MAPPING
      IF(ASSOCIATED(equationsmapping)) THEN
        dynamicmatrices=>equationsmatrices%DYNAMIC_MATRICES
        IF(ASSOCIATED(dynamicmatrices)) THEN
          !Initialise the dynamic nodal matrices
          dynamicmapping=>equationsmapping%DYNAMIC_MAPPING
          IF(ASSOCIATED(dynamicmapping)) THEN
            DO matrixidx=1,dynamicmatrices%NUMBER_OF_DYNAMIC_MATRICES
              equationsmatrix=>dynamicmatrices%MATRICES(matrixidx)%PTR
              IF(ASSOCIATED(equationsmatrix)) THEN
                fieldvariable=>dynamicmapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrixidx)%VARIABLE
                CALL equationsmatrices_nodalmatrixsetup(equationsmatrix%NodalMatrix,fieldvariable,fieldvariable, &
                  & err,error,*999)
              ELSE
                localerror="Equations dynamic matrix number "//trim(numbertovstring(matrixidx,"*",err,error))// &
                  & " is not associated."
                CALL flagerror(localerror,err,error,*999)
              ENDIF
            ENDDO !matrixIdx
          ELSE
            CALL flagerror("Equations mapping dynamic mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        linearmatrices=>equationsmatrices%LINEAR_MATRICES
        IF(ASSOCIATED(linearmatrices)) THEN
          !Initialise the linear nodal matrices
          linearmapping=>equationsmapping%LINEAR_MAPPING
          IF(ASSOCIATED(linearmapping)) THEN
            DO matrixidx=1,linearmatrices%NUMBER_OF_LINEAR_MATRICES
              equationsmatrix=>linearmatrices%MATRICES(matrixidx)%PTR
              IF(ASSOCIATED(equationsmatrix)) THEN
                fieldvariable=>linearmapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrixidx)%VARIABLE
                CALL equationsmatrices_nodalmatrixsetup(equationsmatrix%NodalMatrix,fieldvariable,fieldvariable, &
                  & err,error,*999)
              ELSE
                localerror="Equations linear matrix number "//trim(numbertovstring(matrixidx,"*",err,error))// &
                  & " is not associated."
                CALL flagerror(localerror,err,error,*999)
              ENDIF
            ENDDO !matrixIdx
          ELSE
            CALL flagerror("Equations mapping linear mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        nonlinearmatrices=>equationsmatrices%NONLINEAR_MATRICES
        IF(ASSOCIATED(nonlinearmatrices)) THEN
          !Initialise the Jacobian nodal matrices
          nonlinearmapping=>equationsmapping%NONLINEAR_MAPPING
          IF(ASSOCIATED(nonlinearmapping)) THEN
            fieldvariable=>nonlinearmapping%JACOBIAN_TO_VAR_MAP(1)%VARIABLE
            DO matrixidx=1,nonlinearmatrices%NUMBER_OF_JACOBIANS
              jacobianmatrix=>nonlinearmatrices%JACOBIANS(matrixidx)%PTR
              IF(ASSOCIATED(jacobianmatrix)) THEN
                columnfieldvariable=>nonlinearmapping%JACOBIAN_TO_VAR_MAP(matrixidx)%VARIABLE
                CALL equationsmatrices_nodalmatrixsetup(jacobianmatrix%NodalJacobian,fieldvariable,columnfieldvariable, &
                  & err,error,*999)
              ELSE
                CALL flagerror("Jacobian matrix is not associated.",err,error,*999)
              ENDIF
            ENDDO
            !Use RHS variable for residual vector, otherwise first nonlinear variable if no RHS
            rhsmapping=>equationsmapping%RHS_MAPPING
            IF(ASSOCIATED(rhsmapping)) THEN
              fieldvariable=>rhsmapping%RHS_VARIABLE
            ELSE
              fieldvariable=>nonlinearmapping%JACOBIAN_TO_VAR_MAP(1)%VARIABLE
            ENDIF
            CALL equationsmatrices_nodalvectorsetup(nonlinearmatrices%NodalResidual,fieldvariable,err,error,*999)
            nonlinearmatrices%NodalResidualCalculated=0
          ELSE
            CALL flagerror("Equations mapping nonlinear mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        rhsvector=>equationsmatrices%RHS_VECTOR
        IF(ASSOCIATED(rhsvector)) THEN
          !Initialise the RHS nodal vector
          rhsmapping=>equationsmapping%RHS_MAPPING
          IF(ASSOCIATED(rhsmapping)) THEN
            fieldvariable=>rhsmapping%RHS_VARIABLE
            CALL equationsmatrices_nodalvectorsetup(rhsvector%NodalVector,fieldvariable,err,error,*999)
          ELSE
            CALL flagerror("RHS mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        sourcevector=>equationsmatrices%SOURCE_VECTOR
        IF(ASSOCIATED(sourcevector)) THEN
          !Initialise the source nodal vector. Note that the number of rows in the source vector is taken, for now, from the RHS
          !vector
          IF(ASSOCIATED(rhsvector)) THEN
            !Initialise the RHS nodal vector
            rhsmapping=>equationsmapping%RHS_MAPPING
            IF(ASSOCIATED(rhsmapping)) THEN
              fieldvariable=>rhsmapping%RHS_VARIABLE
              CALL equationsmatrices_nodalvectorsetup(sourcevector%NodalVector,fieldvariable,err,error,*999)
            ELSE
              CALL flagerror("RHS mapping is not associated.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Not implemented.",err,error,*999)
          ENDIF
        ENDIF
      ELSE
        CALL flagerror("Equations matrices mapping is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF
    
    exits("EquationsMatrices_NodalInitialise")
    RETURN
999 errorsexits("EquationsMatrices_NodalInitialise",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_nodalinitialise

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

  SUBROUTINE equationsmatrices_nodalmatrixsetup(nodalMatrix,rowsFieldVariable,colsFieldVariable,err,error,*)

    !Argument variables
    TYPE(nodalmatrixtype) :: nodalMatrix
    TYPE(field_variable_type), POINTER :: rowsFieldVariable
    TYPE(field_variable_type), POINTER :: colsFieldVariable
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: dummyErr
    TYPE(varying_string) :: dummyError

    enters("EquationsMatrices_NodalMatrixSetup",err,error,*998)

    IF(ASSOCIATED(rowsfieldvariable)) THEN
      IF(ASSOCIATED(colsfieldvariable)) THEN
        nodalmatrix%maxNumberOfRows=rowsfieldvariable%maxNumberNodeInterpolationParameters* &
          & rowsfieldvariable%NUMBER_OF_COMPONENTS
        nodalmatrix%maxNumberOfColumns=colsfieldvariable%maxNumberNodeInterpolationParameters* &
          & colsfieldvariable%NUMBER_OF_COMPONENTS
        IF(ALLOCATED(nodalmatrix%rowDofs)) THEN
          CALL flagerror("Nodal matrix row dofs already allocated.",err,error,*999)
        ELSE
          ALLOCATE(nodalmatrix%rowDofs(nodalmatrix%maxNumberOfRows),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate nodal matrix row dofs.",err,error,*999)
        ENDIF
        IF(ALLOCATED(nodalmatrix%columnDofs)) THEN
          CALL flagerror("Nodal matrix column dofs already allocated.",err,error,*999)
        ELSE
          ALLOCATE(nodalmatrix%columnDofs(nodalmatrix%maxNumberOfColumns),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate nodal matrix column dofs.",err,error,*999)
        ENDIF
        IF(ALLOCATED(nodalmatrix%matrix)) THEN
          CALL flagerror("Nodal matrix already allocated.",err,error,*999)
        ELSE
          ALLOCATE(nodalmatrix%matrix(nodalmatrix%maxNumberOfRows,nodalmatrix%maxNumberOfColumns),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate nodal matrix.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Columns field variable is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Rows field variable is not associated.",err,error,*999)
    ENDIF
    
    exits("EquationsMatrices_NodalMatrixSetup")
    RETURN
999 CALL equationsmatrices_nodalmatrixfinalise(nodalmatrix,dummyerr,dummyerror,*998)
998 errorsexits("EquationsMatrices_NodalMatrixSetup",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_nodalmatrixsetup

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

  SUBROUTINE equationsmatrices_nodalvectorsetup(nodalVector,rowsFieldVariable,err,error,*)

    !Argument variables
    TYPE(nodalvectortype) :: nodalVector
    TYPE(field_variable_type), POINTER :: rowsFieldVariable
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: dummyErr
    TYPE(varying_string) :: dummyError

    enters("EquationsMatrices_NodalVectorSetup",err,error,*998)

    IF(ASSOCIATED(rowsfieldvariable)) THEN
      nodalvector%maxNumberOfRows=rowsfieldvariable%maxNumberNodeInterpolationParameters* &
        & rowsfieldvariable%NUMBER_OF_COMPONENTS
      IF(ALLOCATED(nodalvector%rowDofs)) THEN
        CALL flagerror("Nodal vector row dofs is already allocated.",err,error,*999)        
      ELSE
        ALLOCATE(nodalvector%rowDofs(nodalvector%maxNumberOfRows),stat=err)
        IF(err/=0) CALL flagerror("Could not allocate nodal vector row dofs.",err,error,*999)
      ENDIF
      IF(ALLOCATED(nodalvector%vector)) THEN
        CALL flagerror("Nodal vector vector already allocated.",err,error,*999)        
      ELSE
        ALLOCATE(nodalvector%vector(nodalvector%maxNumberOfRows),stat=err)
        IF(err/=0) CALL flagerror("Could not allocate nodal vector vector.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Rows field variable is not associated.",err,error,*999)
    ENDIF
    
    exits("EquationsMatrices_NodalVectorSetup")
    RETURN
999 CALL equationsmatrices_nodalvectorfinalise(nodalvector,dummyerr,dummyerror,*998)
998 errorsexits("EquationsMatrices_NodalVectorSetup",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_nodalvectorsetup

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

  SUBROUTINE equationsmatrices_nodalfinalise(equationsMatrices,err,error,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: equationsMatrices
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: matrixIdx
    TYPE(equations_jacobian_type), POINTER :: jacobianMatrix
    TYPE(equations_matrices_dynamic_type), POINTER :: dynamicMatrices
    TYPE(equations_matrices_linear_type), POINTER :: linearMatrices
    TYPE(equations_matrices_nonlinear_type), POINTER :: nonlinearMatrices
    TYPE(equations_matrices_rhs_type), POINTER :: rhsVector
    TYPE(equations_matrices_source_type), POINTER :: sourceVector
    TYPE(equations_matrix_type), POINTER :: equationsMatrix
    TYPE(varying_string) :: localError
    
    enters("EquationsMatrices_NodalFinalise",err,error,*999)

    IF(ASSOCIATED(equationsmatrices)) THEN
      dynamicmatrices=>equationsmatrices%DYNAMIC_MATRICES
      IF(ASSOCIATED(dynamicmatrices)) THEN
        !Finalise the dynamic nodal matrices
        DO matrixidx=1,dynamicmatrices%NUMBER_OF_DYNAMIC_MATRICES
          equationsmatrix=>dynamicmatrices%MATRICES(matrixidx)%PTR
          IF(ASSOCIATED(equationsmatrix)) THEN
            equationsmatrix%NodalMatrix%maxNumberOfRows=0
            equationsmatrix%NodalMatrix%maxNumberOfColumns=0
            IF(ALLOCATED(equationsmatrix%NodalMatrix%rowDofs)) DEALLOCATE(equationsmatrix%NodalMatrix%rowDofs)
            IF(ALLOCATED(equationsmatrix%NodalMatrix%columnDofs)) DEALLOCATE(equationsmatrix%NodalMatrix%columnDofs)
            IF(ALLOCATED(equationsmatrix%NodalMatrix%matrix)) DEALLOCATE(equationsmatrix%NodalMatrix%matrix)
          ELSE
            localerror="Equations matrix for dynamic matrix number "//trim(numbertovstring(matrixidx,"*",err,error))// &
              & " is not associated."
            CALL flagerror(localerror,err,error,*999)
          ENDIF
        ENDDO !matrixIdx
      ENDIF
      linearmatrices=>equationsmatrices%LINEAR_MATRICES
      IF(ASSOCIATED(linearmatrices)) THEN
        !Finalise the linear nodal matrices
        DO matrixidx=1,linearmatrices%NUMBER_OF_LINEAR_MATRICES
          equationsmatrix=>linearmatrices%MATRICES(matrixidx)%PTR
          IF(ASSOCIATED(equationsmatrix)) THEN
            equationsmatrix%NodalMatrix%maxNumberOfRows=0
            equationsmatrix%NodalMatrix%maxNumberOfColumns=0
            IF(ALLOCATED(equationsmatrix%NodalMatrix%rowDofs)) DEALLOCATE(equationsmatrix%NodalMatrix%rowDofs)
            IF(ALLOCATED(equationsmatrix%NodalMatrix%columnDofs)) DEALLOCATE(equationsmatrix%NodalMatrix%columnDofs)
            IF(ALLOCATED(equationsmatrix%NodalMatrix%matrix)) DEALLOCATE(equationsmatrix%NodalMatrix%matrix)
          ELSE
            localerror="Equations matrix for linear matrix number "//trim(numbertovstring(matrixidx,"*",err,error))// &
              & " is not associated."
            CALL flagerror(localerror,err,error,*999)
          ENDIF
        ENDDO !matrixIdx
      ENDIF
      nonlinearmatrices=>equationsmatrices%NONLINEAR_MATRICES
      IF(ASSOCIATED(nonlinearmatrices)) THEN
        DO matrixidx=1,nonlinearmatrices%NUMBER_OF_JACOBIANS
          jacobianmatrix=>nonlinearmatrices%JACOBIANS(matrixidx)%PTR
          IF(ASSOCIATED(jacobianmatrix)) THEN
            jacobianmatrix%NodalJacobian%maxNumberOfRows=0
            jacobianmatrix%NodalJacobian%maxNumberOfColumns=0
            IF(ALLOCATED(jacobianmatrix%NodalJacobian%rowDofs)) DEALLOCATE(jacobianmatrix%NodalJacobian%rowDofs)
            IF(ALLOCATED(jacobianmatrix%NodalJacobian%columnDofs)) DEALLOCATE(jacobianmatrix%NodalJacobian%columnDofs)
            IF(ALLOCATED(jacobianmatrix%NodalJacobian%matrix)) DEALLOCATE(jacobianmatrix%NodalJacobian%matrix)
          ELSE
            CALL flagerror("Nonlinear matrices Jacobian number "//trim(numbertovstring(matrixidx,"*",err,error))// &
                & " is not associated.",err,error,*999)
          ENDIF
        ENDDO
        nonlinearmatrices%NodalResidual%maxNumberOfRows=0
        IF(ALLOCATED(nonlinearmatrices%NodalResidual%rowDofs)) DEALLOCATE(nonlinearmatrices%NodalResidual%rowDofs)
        IF(ALLOCATED(nonlinearmatrices%NodalResidual%vector)) DEALLOCATE(nonlinearmatrices%NodalResidual%vector)
      ENDIF
      rhsvector=>equationsmatrices%RHS_VECTOR
      IF(ASSOCIATED(rhsvector)) THEN
        !Finalise the nodal vector
        rhsvector%NodalVector%maxNumberOfRows=0
        IF(ALLOCATED(rhsvector%NodalVector%rowDofs)) DEALLOCATE(rhsvector%NodalVector%rowDofs)
        IF(ALLOCATED(rhsvector%NodalVector%vector)) DEALLOCATE(rhsvector%NodalVector%vector)
      ENDIF
      sourcevector=>equationsmatrices%SOURCE_VECTOR
      IF(ASSOCIATED(sourcevector)) THEN
        !Finalise the nodal source vector
        sourcevector%NodalVector%maxNumberOfRows=0
        IF(ALLOCATED(sourcevector%NodalVector%rowDofs)) DEALLOCATE(sourcevector%NodalVector%rowDofs)
        IF(ALLOCATED(sourcevector%NodalVector%vector)) DEALLOCATE(sourcevector%NodalVector%vector)
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF
    
    exits("EquationsMatrices_NodalFinalise")
    RETURN
999 errorsexits("EquationsMatrices_NodalFinalise",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_nodalfinalise

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

  SUBROUTINE equationsmatrices_nodalmatrixinitialise(nodalMatrix,err,error,*)

    !Argument variables
    TYPE(nodalmatrixtype) :: nodalMatrix !The nodal matrix to initialise
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables

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

    nodalmatrix%equationsMatrixNumber=0
    nodalmatrix%numberOfRows=0
    nodalmatrix%numberOfColumns=0
    nodalmatrix%maxNumberOfRows=0
    nodalmatrix%maxNumberOfColumns=0
       
    exits("EquationsMatrices_NodalMatrixInitialise")
    RETURN
999 errorsexits("EquationsMatrices_NodalMatrixInitialise",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_nodalmatrixinitialise

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

  SUBROUTINE equationsmatrices_nodalmatrixfinalise(nodalMatrix,err,error,*)

    !Argument variables
    TYPE(nodalmatrixtype):: nodalMatrix
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    
    enters("EquationsMatrices_NodalMatrixFinalise",err,error,*999)

    IF(ALLOCATED(nodalmatrix%rowDofs)) DEALLOCATE(nodalmatrix%rowDofs)
    IF(ALLOCATED(nodalmatrix%columnDofs)) DEALLOCATE(nodalmatrix%columnDofs)
    IF(ALLOCATED(nodalmatrix%matrix)) DEALLOCATE(nodalmatrix%matrix)
    
    exits("EquationsMatrices_NodalMatrixFinalise")
    RETURN
999 errorsexits("EquationsMatrices_NodalMatrixFinalise",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_nodalmatrixfinalise

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

  SUBROUTINE equationsmatrices_nodalvectorinitialise(nodalVector,err,error,*)

    !Argument variables
    TYPE(nodalvectortype) :: nodalVector !The nodal vector to initialise
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables

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

    nodalvector%numberOfRows=0
    nodalvector%maxNumberOfRows=0
       
    exits("EquationsMatrices_NodalVectorInitialise")
    RETURN
999 errorsexits("EquationsMatrices_NodalVectorInitialise",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_nodalvectorinitialise

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

  SUBROUTINE equationsmatrices_nodalvectorfinalise(nodalVector,err,error,*)

    !Argument variables
    TYPE(nodalvectortype):: nodalVector
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    
    enters("EquationsMatrices_NodalVectorFinalise",err,error,*999)

    IF(ALLOCATED(nodalvector%rowDofs)) DEALLOCATE(nodalvector%rowDofs)
    IF(ALLOCATED(nodalvector%vector)) DEALLOCATE(nodalvector%vector)
    
    exits("EquationsMatrices_NodalVectorFinalise")
    RETURN
999 errorsexits("EquationsMatrices_NodalVectorFinalise",err,error)

    RETURN 1
  END SUBROUTINE equationsmatrices_nodalvectorfinalise

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

  SUBROUTINE equationsmatrices_jacobiannodeadd(equationsMatrices,err,error,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: equationsMatrices
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: jacobianMatrixIdx
    TYPE(equations_jacobian_type), POINTER :: jacobianMatrix
    TYPE(equations_matrices_nonlinear_type), POINTER :: nonlinearMatrices
    TYPE(varying_string) :: localError
    
#ifdef TAUPROF
    CALL tau_static_phase_start("EquationsMatrices_JacobianNodeAdd()")
#endif

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

    IF(ASSOCIATED(equationsmatrices)) THEN
      nonlinearmatrices=>equationsmatrices%NONLINEAR_MATRICES
      IF(ASSOCIATED(nonlinearmatrices)) THEN
        DO jacobianmatrixidx=1,nonlinearmatrices%NUMBER_OF_JACOBIANS
          jacobianmatrix=>nonlinearmatrices%JACOBIANS(jacobianmatrixidx)%PTR
          IF(ASSOCIATED(jacobianmatrix)) THEN
            IF(jacobianmatrix%UPDATE_JACOBIAN) THEN
              !Add in Jacobian element matrices
              CALL distributed_matrix_values_add(jacobianmatrix%jacobian,jacobianmatrix%NodalJacobian%rowDofs(1: &
                & jacobianmatrix%NodalJacobian%numberOfRows),jacobianmatrix%NodalJacobian%columnDofs(1: &
                & jacobianmatrix%NodalJacobian%numberOfColumns),jacobianmatrix%NodalJacobian%matrix(1: &
                & jacobianmatrix%NodalJacobian%numberOfRows,1:jacobianmatrix%NodalJacobian%numberOfColumns), &
                & err,error,*999)
            ENDIF
          ELSE
            localerror="Jacobian matrix for Jacobian matrix index "// &
              & trim(numbertovstring(jacobianmatrixidx,"*",err,error))//" is not associated."
            CALL flagerror(localerror,err,error,*999)
          ENDIF
        ENDDO !jacobianMatrixIdx
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not allocated.",err,error,*999)
    ENDIF
#ifdef TAUPROF
    CALL tau_static_phase_stop("EquationsMatrices_JacobianNodeAdd()")
#endif
    
    exits("EquationsMatrices_JacobianNodeAdd")
    RETURN
999 errorsexits("EquationsMatrices_JacobianNodeAdd",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_jacobiannodeadd

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

  SUBROUTINE equations_matrices_element_finalise(EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
    TYPE(equations_jacobian_type), POINTER :: JACOBIAN_MATRIX
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
    TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
    TYPE(equations_matrices_source_type), POINTER :: SOURCE_VECTOR
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("EQUATIONS_MATRICES_ELEMENT_FINALISE",err,error,*999)

    IF(ASSOCIATED(equations_matrices)) THEN
      dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
      IF(ASSOCIATED(dynamic_matrices)) THEN
        !Finalise the dynamic element matrices
        DO matrix_idx=1,dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES
          equations_matrix=>dynamic_matrices%MATRICES(matrix_idx)%PTR
          IF(ASSOCIATED(equations_matrix)) THEN
            CALL equations_matrices_element_matrix_finalise(equations_matrix%ELEMENT_MATRIX,err,error,*999)
          ELSE
            local_error="Equations matrix for dynamic matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))// &
              & " is not associated."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ENDDO !matrix_idx
      ENDIF
      linear_matrices=>equations_matrices%LINEAR_MATRICES
      IF(ASSOCIATED(linear_matrices)) THEN
        !Finalise the linear element matrices
        DO matrix_idx=1,linear_matrices%NUMBER_OF_LINEAR_MATRICES
          equations_matrix=>linear_matrices%MATRICES(matrix_idx)%PTR
          IF(ASSOCIATED(equations_matrix)) THEN
            CALL equations_matrices_element_matrix_finalise(equations_matrix%ELEMENT_MATRIX,err,error,*999)
          ELSE
            local_error="Equations matrix for linear matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))// &
              & " is not associated."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ENDDO !matrix_idx
      ENDIF
      nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
      IF(ASSOCIATED(nonlinear_matrices)) THEN
        DO matrix_idx=1,nonlinear_matrices%NUMBER_OF_JACOBIANS
          jacobian_matrix=>nonlinear_matrices%JACOBIANS(matrix_idx)%PTR
          IF(ASSOCIATED(jacobian_matrix)) THEN
            jacobian_matrix%ELEMENT_JACOBIAN%MAX_NUMBER_OF_ROWS=0
            jacobian_matrix%ELEMENT_JACOBIAN%MAX_NUMBER_OF_COLUMNS=0
            IF(ALLOCATED(jacobian_matrix%ELEMENT_JACOBIAN%ROW_DOFS)) DEALLOCATE(jacobian_matrix%ELEMENT_JACOBIAN%ROW_DOFS)
            IF(ALLOCATED(jacobian_matrix%ELEMENT_JACOBIAN%COLUMN_DOFS)) DEALLOCATE(jacobian_matrix%ELEMENT_JACOBIAN%COLUMN_DOFS)
            IF(ALLOCATED(jacobian_matrix%ELEMENT_JACOBIAN%MATRIX)) DEALLOCATE(jacobian_matrix%ELEMENT_JACOBIAN%MATRIX)
          ELSE
            CALL flagerror("Nonlinear matrices Jacobian number "//trim(numbertovstring(matrix_idx,"*",err,error))// &
                & " is not associated.",err,error,*999)
          ENDIF
        ENDDO
        nonlinear_matrices%ELEMENT_RESIDUAL%MAX_NUMBER_OF_ROWS=0
        IF(ALLOCATED(nonlinear_matrices%ELEMENT_RESIDUAL%ROW_DOFS)) DEALLOCATE(nonlinear_matrices%ELEMENT_RESIDUAL%ROW_DOFS)
        IF(ALLOCATED(nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR)) DEALLOCATE(nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR)
      ENDIF
      rhs_vector=>equations_matrices%RHS_VECTOR
      IF(ASSOCIATED(rhs_vector)) THEN
        !Finalise the element vector
        rhs_vector%ELEMENT_VECTOR%MAX_NUMBER_OF_ROWS=0
        IF(ALLOCATED(rhs_vector%ELEMENT_VECTOR%ROW_DOFS)) DEALLOCATE(rhs_vector%ELEMENT_VECTOR%ROW_DOFS)
        IF(ALLOCATED(rhs_vector%ELEMENT_VECTOR%VECTOR)) DEALLOCATE(rhs_vector%ELEMENT_VECTOR%VECTOR)
      ENDIF
      source_vector=>equations_matrices%SOURCE_VECTOR
      IF(ASSOCIATED(source_vector)) THEN
        !Finalise the element source vector
        source_vector%ELEMENT_VECTOR%MAX_NUMBER_OF_ROWS=0
        IF(ALLOCATED(source_vector%ELEMENT_VECTOR%ROW_DOFS)) DEALLOCATE(source_vector%ELEMENT_VECTOR%ROW_DOFS)
        IF(ALLOCATED(source_vector%ELEMENT_VECTOR%VECTOR)) DEALLOCATE(source_vector%ELEMENT_VECTOR%VECTOR)
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF
    
    exits("EQUATIONS_MATRICES_ELEMENT_FINALISE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_ELEMENT_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_element_finalise

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

  SUBROUTINE equations_matrices_element_initialise(EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES !The equations matrices to initialise the element information for
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
    INTEGER(INTG) :: rowsNumberOfElements,colsNumberOfElements !Number of elements in the row and col variables whose dofs are present in the element matrix
    TYPE(equations_jacobian_type), POINTER :: JACOBIAN_MATRIX
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_mapping_dynamic_type), POINTER :: DYNAMIC_MAPPING
    TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
    TYPE(equations_mapping_nonlinear_type), POINTER :: NONLINEAR_MAPPING
    TYPE(equations_mapping_rhs_type), POINTER :: RHS_MAPPING
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
    TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
    TYPE(equations_matrices_source_type), POINTER :: SOURCE_VECTOR
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE,COL_FIELD_VARIABLE
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("EQUATIONS_MATRICES_ELEMENT_INITIALISE",err,error,*999)

    IF(ASSOCIATED(equations_matrices)) THEN
      rowsnumberofelements=1
      colsnumberofelements=1
      equations_mapping=>equations_matrices%EQUATIONS_MAPPING
      IF(ASSOCIATED(equations_mapping)) THEN
        dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
        IF(ASSOCIATED(dynamic_matrices)) THEN
          !Initialise the dynamic element matrices
          dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
          IF(ASSOCIATED(dynamic_mapping)) THEN
            DO matrix_idx=1,dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES
              equations_matrix=>dynamic_matrices%MATRICES(matrix_idx)%PTR
              IF(ASSOCIATED(equations_matrix)) THEN
                field_variable=>dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrix_idx)%VARIABLE
                CALL equations_matrices_element_matrix_setup(equations_matrix%ELEMENT_MATRIX,field_variable,field_variable, &
                  & rowsnumberofelements,colsnumberofelements,err,error,*999)
              ELSE
                local_error="Equations dynamic matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))// &
                  & " is not associated."
                CALL flagerror(local_error,err,error,*999)
              ENDIF
            ENDDO !matrix_idx
          ELSE
            CALL flagerror("Equations mapping dynamic mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        linear_matrices=>equations_matrices%LINEAR_MATRICES
        IF(ASSOCIATED(linear_matrices)) THEN
          !Initialise the linear element matrices
          linear_mapping=>equations_mapping%LINEAR_MAPPING
          IF(ASSOCIATED(linear_mapping)) THEN
            DO matrix_idx=1,linear_matrices%NUMBER_OF_LINEAR_MATRICES
              equations_matrix=>linear_matrices%MATRICES(matrix_idx)%PTR
              IF(ASSOCIATED(equations_matrix)) THEN
                field_variable=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrix_idx)%VARIABLE
                CALL equations_matrices_element_matrix_setup(equations_matrix%ELEMENT_MATRIX,field_variable,field_variable, &
                  & rowsnumberofelements,colsnumberofelements,err,error,*999)
              ELSE
                local_error="Equations linear matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))// &
                  & " is not associated."
                CALL flagerror(local_error,err,error,*999)
              ENDIF
            ENDDO !matrix_idx
          ELSE
            CALL flagerror("Equations mapping linear mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
        IF(ASSOCIATED(nonlinear_matrices)) THEN
          !Initialise the Jacobian element matrices
          nonlinear_mapping=>equations_mapping%NONLINEAR_MAPPING
          IF(ASSOCIATED(nonlinear_mapping)) THEN
            field_variable=>nonlinear_mapping%JACOBIAN_TO_VAR_MAP(1)%VARIABLE
            DO matrix_idx=1,nonlinear_matrices%NUMBER_OF_JACOBIANS
              jacobian_matrix=>nonlinear_matrices%JACOBIANS(matrix_idx)%PTR
              IF(ASSOCIATED(jacobian_matrix)) THEN
                col_field_variable=>nonlinear_mapping%JACOBIAN_TO_VAR_MAP(matrix_idx)%VARIABLE
                CALL equations_matrices_element_matrix_setup(jacobian_matrix%ELEMENT_JACOBIAN,field_variable,col_field_variable, &
                  & rowsnumberofelements,colsnumberofelements,err,error,*999)
              ELSE
                CALL flagerror("Jacobian matrix is not associated.",err,error,*999)
              ENDIF
            ENDDO
            !Use RHS variable for residual vector, otherwise first nonlinear variable if no RHS
            rhs_mapping=>equations_mapping%RHS_MAPPING
            IF(ASSOCIATED(rhs_mapping)) THEN
              field_variable=>rhs_mapping%RHS_VARIABLE
            ELSE
              field_variable=>nonlinear_mapping%JACOBIAN_TO_VAR_MAP(1)%VARIABLE
            ENDIF
            CALL equations_matrices_element_vector_setup(nonlinear_matrices%ELEMENT_RESIDUAL,field_variable,err,error,*999)
            nonlinear_matrices%ELEMENT_RESIDUAL_CALCULATED=0
          ELSE
            CALL flagerror("Equations mapping nonlinear mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        rhs_vector=>equations_matrices%RHS_VECTOR
        IF(ASSOCIATED(rhs_vector)) THEN
          !Initialise the RHS element vector
          rhs_mapping=>equations_mapping%RHS_MAPPING
          IF(ASSOCIATED(rhs_mapping)) THEN
            field_variable=>rhs_mapping%RHS_VARIABLE
            CALL equations_matrices_element_vector_setup(rhs_vector%ELEMENT_VECTOR,field_variable,err,error,*999)
          ELSE
            CALL flagerror("RHS mapping is not associated.",err,error,*999)
          ENDIF
        ENDIF
        source_vector=>equations_matrices%SOURCE_VECTOR
        IF(ASSOCIATED(source_vector)) THEN
          !Initialise the source element vector. Note that the number of rows in the source vector is taken, for now, from the RHS
          !vector
          IF(ASSOCIATED(rhs_vector)) THEN
            !Initialise the RHS element vector
            rhs_mapping=>equations_mapping%RHS_MAPPING
            IF(ASSOCIATED(rhs_mapping)) THEN
              field_variable=>rhs_mapping%RHS_VARIABLE
              CALL equations_matrices_element_vector_setup(source_vector%ELEMENT_VECTOR,field_variable,err,error,*999)
            ELSE
              CALL flagerror("RHS mapping is not associated.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Not implemented.",err,error,*999)
          ENDIF
        ENDIF
      ELSE
        CALL flagerror("Equations matrices mapping is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF
    
    exits("EQUATIONS_MATRICES_ELEMENT_INITIALISE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_ELEMENT_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_element_initialise

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

  SUBROUTINE equations_matrix_finalise(EQUATIONS_MATRIX,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    
    enters("EQUATIONS_MATRIX_FINALISE",err,error,*999)

    IF(ASSOCIATED(equations_matrix)) THEN
      IF(ASSOCIATED(equations_matrix%MATRIX)) CALL distributed_matrix_destroy(equations_matrix%MATRIX,err,error,*999)
      CALL equations_matrices_element_matrix_finalise(equations_matrix%ELEMENT_MATRIX,err,error,*999)
      CALL equationsmatrices_nodalmatrixfinalise(equations_matrix%NodalMatrix,err,error,*999)
      IF(ASSOCIATED(equations_matrix%TEMP_VECTOR)) CALL distributed_vector_destroy(equations_matrix%TEMP_VECTOR,err,error,*999)
    ENDIF
    
    exits("EQUATIONS_MATRIX_FINALISE")
    RETURN
999 errorsexits("EQUATIONS_MATRIX_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrix_finalise

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

  SUBROUTINE equations_matrix_dynamic_initialise(DYNAMIC_MATRICES,MATRIX_NUMBER,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    INTEGER(INTG) :: MATRIX_NUMBER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DUMMY_ERR
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_mapping_dynamic_type), POINTER :: DYNAMIC_MAPPING
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    TYPE(varying_string) :: DUMMY_ERROR,LOCAL_ERROR

    enters("EQUATIONS_MATRIX_DYNAMIC_INITIALISE",err,error,*998)

    IF(ASSOCIATED(dynamic_matrices)) THEN
      IF(matrix_number>0.AND.matrix_number<=dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES) THEN
        equations_matrices=>dynamic_matrices%EQUATIONS_MATRICES
        IF(ASSOCIATED(equations_matrices)) THEN
          equations_mapping=>equations_matrices%EQUATIONS_MAPPING
          IF(ASSOCIATED(equations_mapping)) THEN
            dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
            IF(ASSOCIATED(dynamic_mapping)) THEN
              IF(ASSOCIATED(dynamic_matrices%MATRICES(matrix_number)%PTR)) THEN
                local_error="Equations matrix for dynamic matrix number "//trim(numbertovstring(matrix_number,"*",err,error))// &
                & " is already associated."
                CALL flagerror(local_error,err,error,*998)
              ELSE
                ALLOCATE(dynamic_matrices%MATRICES(matrix_number)%PTR,stat=err)
                IF(err/=0) CALL flagerror("Could not allocate equations matrix.",err,error,*999)
                equations_matrix=>dynamic_matrices%MATRICES(matrix_number)%PTR
                equations_matrix%MATRIX_NUMBER=matrix_number
                equations_matrix%DYNAMIC_MATRICES=>dynamic_matrices
                NULLIFY(equations_matrix%LINEAR_MATRICES)
                equations_matrix%STORAGE_TYPE=matrix_block_storage_type
                equations_matrix%STRUCTURE_TYPE=equations_matrix_no_structure
                equations_matrix%LUMPED=.false.
                equations_matrix%UPDATE_MATRIX=.true.
                equations_matrix%FIRST_ASSEMBLY=.true.
                equations_matrix%NUMBER_OF_COLUMNS=dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrix_number)%NUMBER_OF_COLUMNS
                dynamic_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrix_number)%EQUATIONS_MATRIX=>equations_matrix
                NULLIFY(equations_matrix%MATRIX)
                CALL equationsmatrices_elementmatrixinitialise(equations_matrix%ELEMENT_MATRIX,err,error,*999)
                CALL equationsmatrices_nodalmatrixinitialise(equations_matrix%NodalMatrix,err,error,*999)
                NULLIFY(equations_matrix%TEMP_VECTOR)
              ENDIF
            ELSE
              CALL flagerror("Equations mapping dynamic mapping is not associated.",err,error,*998)
            ENDIF
          ELSE
            CALL flagerror("Equations mapping is not associated.",err,error,*998)
          ENDIF
        ELSE
          CALL flagerror("Dynamic matrices equations matrices is not associated.",err,error,*998)
        ENDIF
      ELSE
        local_error="The specified dynamic matrix number of "//trim(numbertovstring(matrix_number,"*",err,error))// &
          & " is invalid. The matrix number must be > 0 and <= "// &
          & trim(numbertovstring(dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES,"*",err,error))//"."
        CALL flagerror(local_error,err,error,*998)
      ENDIF
    ELSE
      CALL flagerror("Dynamic matrices is not associated.",err,error,*998)
    ENDIF
    
    exits("EQUATIONS_MATRIX_DYNAMIC_INITIALISE")
    RETURN
999 CALL equations_matrix_finalise(dynamic_matrices%MATRICES(matrix_number)%PTR,dummy_err,dummy_error,*998)
998 errorsexits("EQUATIONS_MATRIX_DYNAMIC_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrix_dynamic_initialise

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

  SUBROUTINE equations_matrix_linear_initialise(LINEAR_MATRICES,MATRIX_NUMBER,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
    INTEGER(INTG) :: MATRIX_NUMBER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DUMMY_ERR
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    TYPE(varying_string) :: DUMMY_ERROR,LOCAL_ERROR

    enters("EQUATIONS_MATRIX_LINEAR_INITIALISE",err,error,*998)

    IF(ASSOCIATED(linear_matrices)) THEN
      IF(matrix_number>0.AND.matrix_number<=linear_matrices%NUMBER_OF_LINEAR_MATRICES) THEN
        equations_matrices=>linear_matrices%EQUATIONS_MATRICES
        IF(ASSOCIATED(equations_matrices)) THEN
          equations_mapping=>equations_matrices%EQUATIONS_MAPPING
          IF(ASSOCIATED(equations_mapping)) THEN
            linear_mapping=>equations_mapping%LINEAR_MAPPING
            IF(ASSOCIATED(linear_mapping)) THEN
              IF(ASSOCIATED(linear_matrices%MATRICES(matrix_number)%PTR)) THEN
                local_error="Equations matrix for linear matrix number "//trim(numbertovstring(matrix_number,"*",err,error))// &
                & " is already associated."
                CALL flagerror(local_error,err,error,*998)
              ELSE
                ALLOCATE(linear_matrices%MATRICES(matrix_number)%PTR,stat=err)
                IF(err/=0) CALL flagerror("Could not allocate equations matrix.",err,error,*999)
                equations_matrix=>linear_matrices%MATRICES(matrix_number)%PTR
                equations_matrix%MATRIX_NUMBER=matrix_number
                NULLIFY(equations_matrix%DYNAMIC_MATRICES)
                equations_matrix%LINEAR_MATRICES=>linear_matrices
                equations_matrix%STORAGE_TYPE=matrix_block_storage_type
                equations_matrix%STRUCTURE_TYPE=equations_matrix_no_structure
                equations_matrix%LUMPED=.false.
                equations_matrix%UPDATE_MATRIX=.true.
                equations_matrix%FIRST_ASSEMBLY=.true.
                equations_matrix%NUMBER_OF_COLUMNS=linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrix_number)%NUMBER_OF_COLUMNS
                linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrix_number)%EQUATIONS_MATRIX=>equations_matrix
                NULLIFY(equations_matrix%MATRIX)
                CALL equationsmatrices_elementmatrixinitialise(equations_matrix%ELEMENT_MATRIX,err,error,*999)
                CALL equationsmatrices_nodalmatrixinitialise(equations_matrix%NodalMatrix,err,error,*999)
                NULLIFY(equations_matrix%TEMP_VECTOR)
              ENDIF
            ELSE
              CALL flagerror("Equations mapping linear mapping is not associated.",err,error,*998)
            ENDIF
          ELSE
            CALL flagerror("Equations mapping is not associated.",err,error,*998)
          ENDIF
        ELSE
          CALL flagerror("Linear matrices equations matrices is not associated.",err,error,*998)
        ENDIF
      ELSE
        local_error="The specified linear matrix number of "//trim(numbertovstring(matrix_number,"*",err,error))// &
          & " is invalid. The matrix number must be > 0 and <= "// &
          & trim(numbertovstring(linear_matrices%NUMBER_OF_LINEAR_MATRICES,"*",err,error))//"."
        CALL flagerror(local_error,err,error,*998)
      ENDIF
    ELSE
      CALL flagerror("Linear matrices is not associated.",err,error,*998)
    ENDIF
    
    exits("EQUATIONS_MATRIX_LINEAR_INITIALISE")
    RETURN
999 CALL equations_matrix_finalise(linear_matrices%MATRICES(matrix_number)%PTR,dummy_err,dummy_error,*998)
998 errorsexits("EQUATIONS_MATRIX_LINEAR_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrix_linear_initialise

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

  SUBROUTINE equations_matrices_dynamic_finalise(DYNAMIC_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
     
    enters("EQUATIONS_MATRICES_DYNAMIC_FINALISE",err,error,*999)

    IF(ASSOCIATED(dynamic_matrices)) THEN
      IF(ALLOCATED(dynamic_matrices%MATRICES)) THEN
        DO matrix_idx=1,SIZE(dynamic_matrices%MATRICES,1)
          CALL equations_matrix_finalise(dynamic_matrices%MATRICES(matrix_idx)%PTR,err,error,*999)
        ENDDO !matrix_idx
        DEALLOCATE(dynamic_matrices%MATRICES)
      ENDIF
      IF(ASSOCIATED(dynamic_matrices%TEMP_VECTOR)) CALL distributed_vector_destroy(dynamic_matrices%TEMP_VECTOR,err,error,*999)
      DEALLOCATE(dynamic_matrices)
    ENDIF
    
    exits("EQUATIONS_MATRICES_DYNAMIC_FINALISE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_DYNAMIC_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_dynamic_finalise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equations_matrices_dynamic_initialise(EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DUMMY_ERR,matrix_idx
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_mapping_dynamic_type), POINTER :: DYNAMIC_MAPPING
    TYPE(varying_string) :: DUMMY_ERROR
     
    enters("EQUATIONS_MATRICES_DYNAMIC_INITIALISE",err,error,*998)
    
    IF(ASSOCIATED(equations_matrices)) THEN
      equations_mapping=>equations_matrices%EQUATIONS_MAPPING
      IF(ASSOCIATED(equations_mapping)) THEN
        dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
        IF(ASSOCIATED(dynamic_mapping)) THEN
          IF(ASSOCIATED(equations_matrices%DYNAMIC_MATRICES)) THEN
            CALL flagerror("Equations matrices dynamic matrices is already associated.",err,error,*998)
          ELSE
            ALLOCATE(equations_matrices%DYNAMIC_MATRICES,stat=err)
            IF(err/=0) CALL flagerror("Could not allocate equations matrices dynamic matrices.",err,error,*999)
            equations_matrices%DYNAMIC_MATRICES%EQUATIONS_MATRICES=>equations_matrices
            equations_matrices%DYNAMIC_MATRICES%NUMBER_OF_DYNAMIC_MATRICES=dynamic_mapping%NUMBER_OF_DYNAMIC_EQUATIONS_MATRICES
            ALLOCATE(equations_matrices%DYNAMIC_MATRICES%MATRICES(dynamic_mapping%NUMBER_OF_DYNAMIC_EQUATIONS_MATRICES),stat=err)
            IF(err/=0) CALL flagerror("Could not allocate equations matrices dynamic matrices matrices.",err,error,*999)
            DO matrix_idx=1,dynamic_mapping%NUMBER_OF_DYNAMIC_EQUATIONS_MATRICES
              NULLIFY(equations_matrices%DYNAMIC_MATRICES%MATRICES(matrix_idx)%PTR)
              CALL equations_matrix_dynamic_initialise(equations_matrices%DYNAMIC_MATRICES,matrix_idx,err,error,*999)
            ENDDO !matrix_idx
            NULLIFY(equations_matrices%DYNAMIC_MATRICES%TEMP_VECTOR)
          ENDIF
        ENDIF
      ELSE
        CALL flagerror("Equations matrices equations mapping is not associated.",err,error,*998)        
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*998)
    ENDIF
    
    exits("EQUATIONS_MATRICES_DYNAMIC_INITIALISE")
    RETURN
999 CALL equations_matrices_dynamic_finalise(equations_matrices%DYNAMIC_MATRICES,dummy_err,dummy_error,*998)
998 errorsexits("EQUATIONS_MATRICES_DYNAMIC_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_dynamic_initialise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equations_matrices_jacobian_element_add(EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: jacobian_matrix_idx
    TYPE(equations_jacobian_type), POINTER :: JACOBIAN_MATRIX
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
    TYPE(varying_string) :: LOCAL_ERROR
    
#ifdef TAUPROF
    CALL tau_static_phase_start("EQUATIONS_MATRICES_JACOBIAN_ELEMENT_ADD()")
#endif

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

    IF(ASSOCIATED(equations_matrices)) THEN
      nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
      IF(ASSOCIATED(nonlinear_matrices)) THEN
        DO jacobian_matrix_idx=1,nonlinear_matrices%NUMBER_OF_JACOBIANS
          jacobian_matrix=>nonlinear_matrices%JACOBIANS(jacobian_matrix_idx)%PTR
          IF(ASSOCIATED(jacobian_matrix)) THEN
            IF(jacobian_matrix%UPDATE_JACOBIAN) THEN
              !Add in Jacobian element matrices
              CALL distributed_matrix_values_add(jacobian_matrix%JACOBIAN,jacobian_matrix%ELEMENT_JACOBIAN%ROW_DOFS(1: &
                & jacobian_matrix%ELEMENT_JACOBIAN%NUMBER_OF_ROWS),jacobian_matrix%ELEMENT_JACOBIAN%COLUMN_DOFS(1: &
                & jacobian_matrix%ELEMENT_JACOBIAN%NUMBER_OF_COLUMNS),jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(1: &
                & jacobian_matrix%ELEMENT_JACOBIAN%NUMBER_OF_ROWS,1:jacobian_matrix%ELEMENT_JACOBIAN%NUMBER_OF_COLUMNS), &
                & err,error,*999)
            ENDIF
          ELSE
            local_error="Jacobian matrix for Jacobian matrix index "// &
              & trim(numbertovstring(jacobian_matrix_idx,"*",err,error))//" is not associated."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ENDDO !jacobian_matrix_idx
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not allocated.",err,error,*999)
    ENDIF
#ifdef TAUPROF
    CALL tau_static_phase_stop("EQUATIONS_MATRICES_JACOBIAN_ELEMENT_ADD()")
#endif
    
    exits("EQUATIONS_MATRICES_JACOBIAN_ELEMENT_ADD")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_JACOBIAN_ELEMENT_ADD",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_jacobian_element_add

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

  SUBROUTINE equations_matrices_jacobian_output(ID,EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    INTEGER(INTG), INTENT(IN) :: ID
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: jacobian_matrix_idx
    TYPE(equations_jacobian_type), POINTER :: JACOBIAN_MATRIX
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("EQUATIONS_MATRICES_JACOBIAN_OUTPUT",err,error,*999)
    
    IF(ASSOCIATED(equations_matrices)) THEN
      IF(equations_matrices%EQUATIONS_MATRICES_FINISHED) THEN
        nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
        IF(ASSOCIATED(nonlinear_matrices)) THEN
          CALL write_string(id,"",err,error,*999)
          CALL write_string(id,"Jacobian matrices:",err,error,*999)
          DO jacobian_matrix_idx=1,nonlinear_matrices%NUMBER_OF_JACOBIANS
            jacobian_matrix=>nonlinear_matrices%JACOBIANS(jacobian_matrix_idx)%PTR
            IF(ASSOCIATED(jacobian_matrix)) THEN        
              CALL write_string(id,"Jacobian matrix:",err,error,*999)
              CALL distributed_matrix_output(id,jacobian_matrix%JACOBIAN,err,error,*999)
            ELSE
              local_error="Jacobian matrix for Jacobian matrix index "// &
                & trim(numbertovstring(jacobian_matrix_idx,"*",err,error))//" is not associated."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ENDDO !jacobian_matrix_idx
        ENDIF
      ELSE
        CALL flagerror("Equations matrices have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF
    
    exits("EQUATIONS_MATRICES_JACOBIAN_OUTPUT")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_JACOBIAN_OUTPUT",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_jacobian_output
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equationsmatrices_jacobiantypesset(equationsMatrices,jacobianTypes,err,error,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: equationsMatrices
    INTEGER(INTG), INTENT(IN) :: jacobianTypes(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(equations_matrices_nonlinear_type), POINTER :: nonlinearMatrices
    INTEGER(INTG) :: matrixIdx,numberOfjacobians,jacobianType
    TYPE(varying_string) :: localError

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

    IF(ASSOCIATED(equationsmatrices)) THEN
      nonlinearmatrices=>equationsmatrices%NONLINEAR_MATRICES
      IF(ASSOCIATED(nonlinearmatrices)) THEN
        numberofjacobians=SIZE(jacobiantypes,1)
        IF(numberofjacobians==nonlinearmatrices%NUMBER_OF_JACOBIANS) THEN
          DO matrixidx=1,numberofjacobians
            jacobiantype=jacobiantypes(matrixidx)
            SELECT CASE(jacobiantype)
            CASE(equations_jacobian_finite_difference_calculated, &
                & equations_jacobian_analytic_calculated)
              nonlinearmatrices%JACOBIANS(matrixidx)%PTR%JACOBIAN_CALCULATION_TYPE=jacobiantype
            CASE DEFAULT
              localerror="Invalid Jacobian calculation type of " &
                & //trim(numbertovstring(jacobiantype,"*",err,error))//"."
              CALL flagerror(localerror,err,error,*999)
            END SELECT
          END DO
        ELSE
          localerror="Invalid number of Jacobian calculation types. The number of types " &
            & //trim(numbertovstring(numberofjacobians,"*",err,error)) &
            & //" should be "//trim(numbertovstring(nonlinearmatrices%NUMBER_OF_JACOBIANS,"*",err,error))
          CALL flagerror(localerror,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Equations matrices nonlinear matrices are not associated",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations matrices are not associated",err,error,*999)
    ENDIF

    exits("EquationsMatrices_JacobianTypesSet")
    RETURN
999 errorsexits("EquationsMatrices_JacobianTypesSet",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_jacobiantypesset

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

  SUBROUTINE equations_matrices_linear_finalise(LINEAR_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
     
    enters("EQUATIONS_MATRICES_LINEAR_FINALISE",err,error,*999)

    IF(ASSOCIATED(linear_matrices)) THEN
      IF(ALLOCATED(linear_matrices%MATRICES)) THEN
        DO matrix_idx=1,SIZE(linear_matrices%MATRICES,1)
          CALL equations_matrix_finalise(linear_matrices%MATRICES(matrix_idx)%PTR,err,error,*999)
        ENDDO !matrix_idx
        DEALLOCATE(linear_matrices%MATRICES)
      ENDIF
      DEALLOCATE(linear_matrices)
    ENDIF
    
    exits("EQUATIONS_MATRICES_LINEAR_FINALISE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_LINEAR_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_linear_finalise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equations_matrices_linear_initialise(EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DUMMY_ERR,matrix_idx
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
    TYPE(varying_string) :: DUMMY_ERROR
     
    enters("EQUATIONS_MATRICES_LINEAR_INITIALISE",err,error,*998)
    
    IF(ASSOCIATED(equations_matrices)) THEN
      equations_mapping=>equations_matrices%EQUATIONS_MAPPING
      IF(ASSOCIATED(equations_mapping)) THEN
        linear_mapping=>equations_mapping%LINEAR_MAPPING
        IF(ASSOCIATED(linear_mapping)) THEN
          IF(ASSOCIATED(equations_matrices%LINEAR_MATRICES)) THEN
            CALL flagerror("Equations matrices linear matrices is already associated.",err,error,*998)
          ELSE
            ALLOCATE(equations_matrices%LINEAR_MATRICES,stat=err)
            IF(err/=0) CALL flagerror("Could not allocate equations matrices linear matrices.",err,error,*999)
            equations_matrices%LINEAR_MATRICES%EQUATIONS_MATRICES=>equations_matrices
            equations_matrices%LINEAR_MATRICES%NUMBER_OF_LINEAR_MATRICES=linear_mapping%NUMBER_OF_LINEAR_EQUATIONS_MATRICES
            ALLOCATE(equations_matrices%LINEAR_MATRICES%MATRICES(linear_mapping%NUMBER_OF_LINEAR_EQUATIONS_MATRICES),stat=err)
            IF(err/=0) CALL flagerror("Could not allocate equations matrices linear matrices matrices.",err,error,*999)
            DO matrix_idx=1,linear_mapping%NUMBER_OF_LINEAR_EQUATIONS_MATRICES
              NULLIFY(equations_matrices%LINEAR_MATRICES%MATRICES(matrix_idx)%PTR)
              CALL equations_matrix_linear_initialise(equations_matrices%LINEAR_MATRICES,matrix_idx,err,error,*999)
            ENDDO !matrix_idx
          ENDIF
        ENDIF
      ELSE
        CALL flagerror("Equations matrices equations mapping is not associated.",err,error,*998)        
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*998)
    ENDIF
    
    exits("EQUATIONS_MATRICES_LINEAR_INITIALISE")
    RETURN
999 CALL equations_matrices_linear_finalise(equations_matrices%LINEAR_MATRICES,dummy_err,dummy_error,*998)
998 errorsexits("EQUATIONS_MATRICES_LINEAR_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_linear_initialise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equations_matrices_nonlinear_finalise(NONLINEAR_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
     
    enters("EQUATIONS_MATRICES_NONLINEAR_FINALISE",err,error,*999)

    IF(ASSOCIATED(nonlinear_matrices)) THEN
      IF(ALLOCATED(nonlinear_matrices%JACOBIANS)) THEN
        DO matrix_idx=1,nonlinear_matrices%NUMBER_OF_JACOBIANS
          CALL equations_jacobian_finalise(nonlinear_matrices%JACOBIANS(matrix_idx)%PTR,err,error,*999)
        ENDDO
        DEALLOCATE(nonlinear_matrices%JACOBIANS)
      ENDIF
      IF(ASSOCIATED(nonlinear_matrices%RESIDUAL)) CALL distributed_vector_destroy(nonlinear_matrices%RESIDUAL,err,error,*999)
      CALL equations_matrices_element_vector_finalise(nonlinear_matrices%ELEMENT_RESIDUAL,err,error,*999)
      CALL equationsmatrices_nodalvectorfinalise(nonlinear_matrices%NodalResidual,err,error,*999)
      DEALLOCATE(nonlinear_matrices)
    ENDIF
    
    exits("EQUATIONS_MATRICES_NONLINEAR_FINALISE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_NONLINEAR_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_nonlinear_finalise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equations_matrices_nonlinear_initialise(EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx,DUMMY_ERR
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_mapping_nonlinear_type), POINTER :: NONLINEAR_MAPPING
    TYPE(varying_string) :: DUMMY_ERROR
    
    enters("EQUATIONS_MATRICES_NONLINEAR_INITIALISE",err,error,*998)

    IF(ASSOCIATED(equations_matrices)) THEN
      equations_mapping=>equations_matrices%EQUATIONS_MAPPING
      IF(ASSOCIATED(equations_mapping)) THEN
        nonlinear_mapping=>equations_mapping%NONLINEAR_MAPPING
        IF(ASSOCIATED(nonlinear_mapping)) THEN
          IF(ASSOCIATED(equations_matrices%NONLINEAR_MATRICES)) THEN
            CALL flagerror("Equations matrices nonlinear matrices is already associated.",err,error,*998)
          ELSE
            ALLOCATE(equations_matrices%NONLINEAR_MATRICES,stat=err)
            IF(err/=0) CALL flagerror("Could not allocate equations matrices nonlinear matrices.",err,error,*999)
            equations_matrices%NONLINEAR_MATRICES%EQUATIONS_MATRICES=>equations_matrices
            equations_matrices%NONLINEAR_MATRICES%UPDATE_RESIDUAL=.true.
            equations_matrices%NONLINEAR_MATRICES%FIRST_ASSEMBLY=.true.
            NULLIFY(equations_matrices%NONLINEAR_MATRICES%RESIDUAL)
            CALL equationsmatrices_elementvectorinitialise(equations_matrices%NONLINEAR_MATRICES%ELEMENT_RESIDUAL,err,error,*999)
            CALL equationsmatrices_nodalvectorinitialise(equations_matrices%NONLINEAR_MATRICES%NodalResidual,err,error,*999)
            equations_matrices%NONLINEAR_MATRICES%NUMBER_OF_JACOBIANS=nonlinear_mapping%NUMBER_OF_RESIDUAL_VARIABLES
            ALLOCATE(equations_matrices%NONLINEAR_MATRICES%JACOBIANS(nonlinear_mapping%NUMBER_OF_RESIDUAL_VARIABLES),stat=err)
            IF(err/=0) CALL flagerror("Could not allocate equations matrices Jacobian matrices.",err,error,*999)
            DO matrix_idx=1,nonlinear_mapping%NUMBER_OF_RESIDUAL_VARIABLES
              NULLIFY(equations_matrices%NONLINEAR_MATRICES%JACOBIANS(matrix_idx)%PTR)
              CALL equations_jacobian_initialise(equations_matrices%NONLINEAR_MATRICES,matrix_idx,err,error,*999)
            ENDDO !matrix_idx
          ENDIF
        ENDIF
      ELSE
        CALL flagerror("Equations matrices equations mapping is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*998)
    ENDIF
    
    exits("EQUATIONS_MATRICES_NONLINEAR_INITIALISE")
    RETURN
999 CALL equations_matrices_nonlinear_finalise(equations_matrices%NONLINEAR_MATRICES,dummy_err,dummy_error,*998)
998 errorsexits("EQUATIONS_MATRICES_NONLINEAR_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_nonlinear_initialise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equations_matrices_output(ID,EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    INTEGER(INTG), INTENT(IN) :: ID
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
    TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
    TYPE(equations_matrices_source_type), POINTER :: SOURCE_VECTOR
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    
    enters("EQUATIONS_MATRICES_OUTPUT",err,error,*999)
    
    IF(ASSOCIATED(equations_matrices)) THEN
      IF(equations_matrices%EQUATIONS_MATRICES_FINISHED) THEN
        CALL write_string(id,"",err,error,*999)
        CALL write_string(id,"Equations matrices:",err,error,*999)
        dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
        IF(ASSOCIATED(dynamic_matrices)) THEN
          CALL write_string(id,"Dynamic matrices:",err,error,*999)
          CALL write_string_value(id,"Number of dynamic matrices = ",dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES,err,error,*999)
          DO matrix_idx=1,dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES
            equations_matrix=>dynamic_matrices%MATRICES(matrix_idx)%PTR
            IF(ASSOCIATED(equations_matrix)) THEN
              CALL write_string_value(id,"Equations matrix : ",matrix_idx,err,error,*999)
              CALL distributed_matrix_output(id,equations_matrix%MATRIX,err,error,*999)
            ELSE
              CALL flagerror("Equations matrix is not associated.",err,error,*999)
            ENDIF
          ENDDO !matrix_idx
        ENDIF
        linear_matrices=>equations_matrices%LINEAR_MATRICES
        IF(ASSOCIATED(linear_matrices)) THEN
          CALL write_string(id,"Linear matrices:",err,error,*999)
          CALL write_string_value(id,"Number of linear matrices = ",linear_matrices%NUMBER_OF_LINEAR_MATRICES,err,error,*999)
          DO matrix_idx=1,linear_matrices%NUMBER_OF_LINEAR_MATRICES
            equations_matrix=>linear_matrices%MATRICES(matrix_idx)%PTR
            IF(ASSOCIATED(equations_matrix)) THEN
              CALL write_string_value(id,"Equations matrix : ",matrix_idx,err,error,*999)
              CALL distributed_matrix_output(id,equations_matrix%MATRIX,err,error,*999)
            ELSE
              CALL flagerror("Equations matrix is not associated.",err,error,*999)
            ENDIF
          ENDDO !matrix_idx
        ENDIF
        nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
        IF(ASSOCIATED(nonlinear_matrices)) THEN
          CALL write_string(id,"Nonlinear vectors:",err,error,*999)
          IF(ASSOCIATED(nonlinear_matrices%RESIDUAL)) THEN
            CALL write_string(id,"Residual vector:",err,error,*999)
            CALL distributed_vector_output(id,nonlinear_matrices%RESIDUAL,err,error,*999)
          ELSE
            CALL flagerror("Nonlinear matrices residual is not associated.",err,error,*999)
          ENDIF
        ENDIF
        rhs_vector=>equations_matrices%RHS_VECTOR
        IF(ASSOCIATED(rhs_vector)) THEN
          CALL write_string(id,"RHS vector:",err,error,*999)
          CALL distributed_vector_output(id,rhs_vector%VECTOR,err,error,*999)
        ENDIF
        source_vector=>equations_matrices%SOURCE_VECTOR
        IF(ASSOCIATED(source_vector)) THEN
          CALL write_string(id,"Source vector:",err,error,*999)
          CALL distributed_vector_output(id,source_vector%VECTOR,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Equations matrices have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF
    
    exits("EQUATIONS_MATRICES_OUTPUT")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_OUTPUT",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_output
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equations_matrices_rhs_finalise(RHS_VECTOR,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
     
    enters("EQUATIONS_MATRICES_RHS_FINALISE",err,error,*999)

    IF(ASSOCIATED(rhs_vector)) THEN
      IF(ASSOCIATED(rhs_vector%VECTOR)) CALL distributed_vector_destroy(rhs_vector%VECTOR,err,error,*999)
      CALL equations_matrices_element_vector_finalise(rhs_vector%ELEMENT_VECTOR,err,error,*999)
      CALL equationsmatrices_nodalvectorfinalise(rhs_vector%NodalVector,err,error,*999)
      DEALLOCATE(rhs_vector)
    ENDIF      
     
    exits("EQUATIONS_MATRICES_RHS_FINALISE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_RHS_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_rhs_finalise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equations_matrices_rhs_initialise(EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DUMMY_ERR
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_mapping_rhs_type), POINTER :: RHS_MAPPING
    TYPE(varying_string) :: DUMMY_ERROR
    
    enters("EQUATIONS_MATRICES_RHS_INITIALISE",err,error,*998)

    IF(ASSOCIATED(equations_matrices)) THEN
      equations_mapping=>equations_matrices%EQUATIONS_MAPPING
      IF(ASSOCIATED(equations_mapping)) THEN
        rhs_mapping=>equations_mapping%RHS_MAPPING
        IF(ASSOCIATED(rhs_mapping)) THEN
          IF(ASSOCIATED(equations_matrices%RHS_VECTOR)) THEN
            CALL flagerror("Equations matrices RHS vector is already associated.",err,error,*998)
          ELSE
            ALLOCATE(equations_matrices%RHS_VECTOR,stat=err)
            IF(err/=0) CALL flagerror("Could not allocate equations matrices RHS vector.",err,error,*999)
            equations_matrices%RHS_VECTOR%UPDATE_VECTOR=.true.
            equations_matrices%RHS_VECTOR%FIRST_ASSEMBLY=.true.
            NULLIFY(equations_matrices%RHS_VECTOR%VECTOR)
            CALL equationsmatrices_elementvectorinitialise(equations_matrices%RHS_VECTOR%ELEMENT_VECTOR,err,error,*999)
            CALL equationsmatrices_nodalvectorinitialise(equations_matrices%RHS_VECTOR%NodalVector,err,error,*999)
          ENDIF
        ENDIF
      ELSE
        CALL flagerror("Equations matrices equation mapping is not associated.",err,error,*998)
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*998)
    ENDIF
    
    exits("EQUATIONS_MATRICES_RHS_INITIALISE")
    RETURN
999 CALL equations_matrices_rhs_finalise(equations_matrices%RHS_VECTOR,dummy_err,dummy_error,*998)
998 errorsexits("EQUATIONS_MATRICES_RHS_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_rhs_initialise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equations_matrices_source_finalise(SOURCE_VECTOR,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_source_type), POINTER :: SOURCE_VECTOR
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
     
    enters("EQUATIONS_MATRICES_SOURCE_FINALISE",err,error,*999)

    IF(ASSOCIATED(source_vector)) THEN
      IF(ASSOCIATED(source_vector%VECTOR)) CALL distributed_vector_destroy(source_vector%VECTOR,err,error,*999)
      CALL equations_matrices_element_vector_finalise(source_vector%ELEMENT_VECTOR,err,error,*999)
      CALL equationsmatrices_nodalvectorfinalise(source_vector%NodalVector,err,error,*999)
      DEALLOCATE(source_vector)
    ENDIF      
     
    exits("EQUATIONS_MATRICES_SOURCE_FINALISE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_SOURCE_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_source_finalise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equations_matrices_source_initialise(EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DUMMY_ERR
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_mapping_source_type), POINTER :: SOURCE_MAPPING
    TYPE(varying_string) :: DUMMY_ERROR
    
    enters("EQUATIONS_MATRICES_SOURCE_INITIALISE",err,error,*998)

    IF(ASSOCIATED(equations_matrices)) THEN
      equations_mapping=>equations_matrices%EQUATIONS_MAPPING
      IF(ASSOCIATED(equations_mapping)) THEN
        source_mapping=>equations_mapping%SOURCE_MAPPING
        IF(ASSOCIATED(source_mapping)) THEN
          IF(ASSOCIATED(equations_matrices%SOURCE_VECTOR)) THEN
            CALL flagerror("Equations matrices source vector is already associated.",err,error,*998)
          ELSE
            ALLOCATE(equations_matrices%SOURCE_VECTOR,stat=err)
            IF(err/=0) CALL flagerror("Could not allocate equations matrices source vector.",err,error,*999)
            equations_matrices%SOURCE_VECTOR%UPDATE_VECTOR=.true.
            equations_matrices%SOURCE_VECTOR%FIRST_ASSEMBLY=.true.
            NULLIFY(equations_matrices%SOURCE_VECTOR%VECTOR)
            CALL equationsmatrices_elementvectorinitialise(equations_matrices%SOURCE_VECTOR%ELEMENT_VECTOR,err,error,*999)
            CALL equationsmatrices_nodalvectorinitialise(equations_matrices%SOURCE_VECTOR%NodalVector,err,error,*999)
          ENDIF
        ENDIF
      ELSE
        CALL flagerror("Equations matrices equation mapping is not associated.",err,error,*998)
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*998)
    ENDIF
    
    exits("EQUATIONS_MATRICES_SOURCE_INITIALISE")
    RETURN
999 CALL equations_matrices_source_finalise(equations_matrices%SOURCE_VECTOR,dummy_err,dummy_error,*998)
998 errorsexits("EQUATIONS_MATRICES_SOURCE_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_source_initialise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equations_matrices_dynamic_lumping_type_set(EQUATIONS_MATRICES,LUMPING_TYPE,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(IN) :: LUMPING_TYPE(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("EQUATIONS_MATRICES_DYNAMIC_LUMPING_TYPE_SET",err,error,*999)

    IF(ASSOCIATED(equations_matrices)) THEN
      IF(equations_matrices%EQUATIONS_MATRICES_FINISHED) THEN
        CALL flagerror("Equations matrices have already been finished.",err,error,*999)
      ELSE
        dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
        IF(ASSOCIATED(dynamic_matrices)) THEN
          IF(SIZE(lumping_type,1)==dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES) THEN
            DO matrix_idx=1,dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES
              equations_matrix=>dynamic_matrices%MATRICES(matrix_idx)%PTR
              IF(ASSOCIATED(equations_matrix)) THEN
                SELECT CASE(lumping_type(matrix_idx))
                CASE(equations_matrix_unlumped)
                  equations_matrix%LUMPED=.false.
                CASE(equations_matrix_lumped)
                  equations_matrix%LUMPED=.true.        
                CASE DEFAULT
                  local_error="The specified lumping type of "//trim(numbertovstring(lumping_type(matrix_idx),"*",err,error))// &
                    & " for the dynamic matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))//" is invalid."
                  CALL flagerror(local_error,err,error,*999)
                END SELECT
              ELSE
                CALL flagerror("Equations matrix is not associated.",err,error,*999)
              ENDIF
            ENDDO !matrix_idx
          ELSE
            local_error="The size of the lumping type array ("//trim(numbertovstring(SIZE(lumping_type,1),"*",err,error))// &
              & ") is not equal to the number of dynamic matrices ("// &
              & trim(numbertovstring(dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES,"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Equations matrices dynamic matrices is not associated.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF
    
    exits("EQUATIONS_MATRICES_DYNAMIC_LUMPING_TYPE_SET")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_DYNAMIC_LUMPING_TYPE_SET",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_dynamic_lumping_type_set

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

  SUBROUTINE equations_matrices_dynamic_storage_type_set(EQUATIONS_MATRICES,STORAGE_TYPE,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(IN) :: STORAGE_TYPE(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("EQUATIONS_MATRICES_DYNAMIC_STORAGE_TYPE_SET",err,error,*999)

    IF(ASSOCIATED(equations_matrices)) THEN
      IF(equations_matrices%EQUATIONS_MATRICES_FINISHED) THEN
        CALL flagerror("Equations matrices have already been finished.",err,error,*999)
      ELSE
        dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
        IF(ASSOCIATED(dynamic_matrices)) THEN
          IF(SIZE(storage_type,1)==dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES) THEN
            DO matrix_idx=1,dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES
              equations_matrix=>dynamic_matrices%MATRICES(matrix_idx)%PTR
              IF(ASSOCIATED(equations_matrix)) THEN
                SELECT CASE(storage_type(matrix_idx))
                CASE(distributed_matrix_block_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_block_storage_type
                CASE(distributed_matrix_diagonal_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_diagonal_storage_type        
                CASE(distributed_matrix_column_major_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_column_major_storage_type
                CASE(distributed_matrix_row_major_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_row_major_storage_type
                CASE(distributed_matrix_compressed_row_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_compressed_row_storage_type
                CASE(distributed_matrix_compressed_column_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_compressed_column_storage_type
                CASE(distributed_matrix_row_column_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_row_column_storage_type
                CASE DEFAULT
                  local_error="The specified storage type of "//trim(numbertovstring(storage_type(matrix_idx),"*",err,error))// &
                    & " for the dynamic matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))//" is invalid."
                  CALL flagerror(local_error,err,error,*999)
                END SELECT
              ELSE
                CALL flagerror("Equations matrix is not associated.",err,error,*999)
              ENDIF
            ENDDO !matrix_idx
          ELSE
            local_error="The size of the storage type array ("//trim(numbertovstring(SIZE(storage_type,1),"*",err,error))// &
              & ") is not equal to the number of dynamic matrices ("// &
              & trim(numbertovstring(dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES,"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Equations matrices dynamic matrices is not associated.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF
    
    exits("EQUATIONS_MATRICES_DYNAMIC_STORAGE_TYPE_SET")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_DYNAMIC_STORAGE_TYPE_SET",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_dynamic_storage_type_set

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

  SUBROUTINE equations_matrices_linear_storage_type_set(EQUATIONS_MATRICES,STORAGE_TYPE,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(IN) :: STORAGE_TYPE(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
    TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("EQUATIONS_MATRICES_LINEAR_STORAGE_TYPE_SET",err,error,*999)

    IF(ASSOCIATED(equations_matrices)) THEN
      IF(equations_matrices%EQUATIONS_MATRICES_FINISHED) THEN
        CALL flagerror("Equations matrices have been finished.",err,error,*999)
      ELSE
        linear_matrices=>equations_matrices%LINEAR_MATRICES
        IF(ASSOCIATED(linear_matrices)) THEN
          IF(SIZE(storage_type,1)==linear_matrices%NUMBER_OF_LINEAR_MATRICES) THEN
            DO matrix_idx=1,linear_matrices%NUMBER_OF_LINEAR_MATRICES
              equations_matrix=>linear_matrices%MATRICES(matrix_idx)%PTR
              IF(ASSOCIATED(equations_matrix)) THEN
                SELECT CASE(storage_type(matrix_idx))
                CASE(distributed_matrix_block_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_block_storage_type
                CASE(distributed_matrix_diagonal_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_diagonal_storage_type        
                CASE(distributed_matrix_column_major_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_column_major_storage_type
                CASE(distributed_matrix_row_major_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_row_major_storage_type
                CASE(distributed_matrix_compressed_row_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_compressed_row_storage_type
                CASE(distributed_matrix_compressed_column_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_compressed_column_storage_type
                CASE(distributed_matrix_row_column_storage_type)
                  equations_matrix%STORAGE_TYPE=distributed_matrix_row_column_storage_type
                CASE DEFAULT
                  local_error="The specified storage type of "//trim(numbertovstring(storage_type(matrix_idx),"*",err,error))// &
                    & " for the linear matrix number "//trim(numbertovstring(matrix_idx,"*",err,error))//" is invalid."
                  CALL flagerror(local_error,err,error,*999)
                END SELECT
              ELSE
                CALL flagerror("Equations matrix is not associated.",err,error,*999)
              ENDIF
            ENDDO !matrix_idx
          ELSE
            local_error="The size of the storage type array ("//trim(numbertovstring(SIZE(storage_type,1),"*",err,error))// &
              & ") is not equal to the number of linear matrices ("// &
              & trim(numbertovstring(linear_matrices%NUMBER_OF_LINEAR_MATRICES,"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Equations matrices linear matrices is not associated.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF
    
    exits("EQUATIONS_MATRICES_LINEAR_STORAGE_TYPE_SET")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_LINEAR_STORAGE_TYPE_SET",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_linear_storage_type_set

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

  SUBROUTINE equationsmatrices_nonlinearstoragetypeset0(EQUATIONS_MATRICES,STORAGE_TYPE,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(IN) :: STORAGE_TYPE(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
    TYPE(equations_jacobian_type), POINTER :: JACOBIAN_MATRIX
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("EquationsMatrices_NonlinearStorageTypeSet0",err,error,*999)

    IF(ASSOCIATED(equations_matrices)) THEN
      IF(equations_matrices%EQUATIONS_MATRICES_FINISHED) THEN
        CALL flagerror("Equations matrices have been finished.",err,error,*999)
      ELSE
        nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
        IF(ASSOCIATED(nonlinear_matrices)) THEN
          IF(SIZE(storage_type,1)==nonlinear_matrices%NUMBER_OF_JACOBIANS) THEN
            DO matrix_idx=1,nonlinear_matrices%NUMBER_OF_JACOBIANS
              jacobian_matrix=>nonlinear_matrices%JACOBIANS(matrix_idx)%PTR
              IF(ASSOCIATED(jacobian_matrix)) THEN
                SELECT CASE(storage_type(matrix_idx))
                CASE(distributed_matrix_block_storage_type)
                  jacobian_matrix%STORAGE_TYPE=distributed_matrix_block_storage_type
                CASE(distributed_matrix_diagonal_storage_type)
                  jacobian_matrix%STORAGE_TYPE=distributed_matrix_diagonal_storage_type        
                CASE(distributed_matrix_column_major_storage_type)
                  jacobian_matrix%STORAGE_TYPE=distributed_matrix_column_major_storage_type
                CASE(distributed_matrix_row_major_storage_type)
                  jacobian_matrix%STORAGE_TYPE=distributed_matrix_row_major_storage_type
                CASE(distributed_matrix_compressed_row_storage_type)
                  jacobian_matrix%STORAGE_TYPE=distributed_matrix_compressed_row_storage_type
                CASE(distributed_matrix_compressed_column_storage_type)
                  jacobian_matrix%STORAGE_TYPE=distributed_matrix_compressed_column_storage_type
                CASE(distributed_matrix_row_column_storage_type)
                  jacobian_matrix%STORAGE_TYPE=distributed_matrix_row_column_storage_type
                CASE DEFAULT
                  local_error="The specified storage type of "//trim(numbertovstring(storage_type(matrix_idx),"*",err,error))// &
                    & " for the Jacobian matrix is invalid."
                  CALL flagerror(local_error,err,error,*999)
                END SELECT
              ELSE
                CALL flagerror("Jacobian matrix is not associated.",err,error,*999)
              ENDIF
            ENDDO
          ELSE
            local_error="The size of the storage type array ("//trim(numbertovstring(SIZE(storage_type,1),"*",err,error))// &
              & ") is not equal to the number of Jacobian matrices ("// &
              & trim(numbertovstring(nonlinear_matrices%NUMBER_OF_JACOBIANS,"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Equations matrices nonlinear matrices is not associated.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF
    
    exits("EquationsMatrices_NonlinearStorageTypeSet0")
    RETURN
999 errorsexits("EquationsMatrices_NonlinearStorageTypeSet0",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrices_nonlinearstoragetypeset0

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

  SUBROUTINE equationsmatrices_nonlinearstoragetypeset1(EQUATIONS_MATRICES,STORAGE_TYPE,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(IN) :: STORAGE_TYPE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG), ALLOCATABLE :: STORAGE_TYPES(:)
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES

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

    IF(ASSOCIATED(equations_matrices)) THEN
      IF(equations_matrices%EQUATIONS_MATRICES_FINISHED) THEN
        CALL flagerror("Equations matrices have been finished.",err,error,*999)
      ELSE
        nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
        IF(ASSOCIATED(nonlinear_matrices)) THEN
          ALLOCATE(storage_types(nonlinear_matrices%NUMBER_OF_JACOBIANS),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate storage types.",err,error,*999)
          storage_types=storage_type
          CALL equationsmatrices_nonlinearstoragetypeset0(equations_matrices,storage_types,err,error,*999)
          DEALLOCATE(storage_types)
        ELSE
          CALL flagerror("Equations matrices nonlinear matrices is not associated.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF

    exits("EquationsMatrices_NonlinearStorageTypeSet1")
    RETURN
999 errorsexits("EquationsMatrices_NonlinearStorageTypeSet1",err,error)
    RETURN 1
    
  END SUBROUTINE equationsmatrices_nonlinearstoragetypeset1

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

  SUBROUTINE equationsmatrices_dynamicstructuretypeset(EQUATIONS_MATRICES,STRUCTURE_TYPE,ERR,ERROR,*)
    
    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(IN) :: STRUCTURE_TYPE(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    TYPE(varying_string) :: LOCAL_ERROR

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

    IF(ASSOCIATED(equations_matrices)) THEN
      IF(equations_matrices%EQUATIONS_MATRICES_FINISHED) THEN
        CALL flagerror("Equations matrices have been finished.",err,error,*999)
      ELSE
        dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
        IF(ASSOCIATED(dynamic_matrices)) THEN
          IF(SIZE(structure_type,1)==dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES) THEN
           DO matrix_idx=1,dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES
              equations_matrix=>dynamic_matrices%MATRICES(matrix_idx)%PTR
              IF(ASSOCIATED(equations_matrix)) THEN
                SELECT CASE(structure_type(matrix_idx))
                CASE(equations_matrix_no_structure)
                  equations_matrix%STRUCTURE_TYPE=equations_matrix_no_structure
                CASE(equations_matrix_fem_structure)
                  equations_matrix%STRUCTURE_TYPE=equations_matrix_fem_structure
                CASE(equations_matrix_diagonal_structure)
                  equations_matrix%STRUCTURE_TYPE=equations_matrix_diagonal_structure
                CASE(equations_matrix_nodal_structure)
                  equations_matrix%STRUCTURE_TYPE=equations_matrix_nodal_structure
                CASE DEFAULT
                  local_error="The specified strucutre type of "// &
                    & trim(numbertovstring(structure_type(matrix_idx),"*",err,error))//" for dynamic matrix number "// &
                    & trim(numbertovstring(matrix_idx,"*",err,error))//" is invalid."
                  CALL flagerror(local_error,err,error,*999)
                END SELECT
              ELSE
                CALL flagerror("Equations matrix is not associated.",err,error,*999)
              ENDIF
            ENDDO !matrix_idx
          ELSE
            local_error="The size of the structure type array ("//trim(numbertovstring(SIZE(structure_type,1),"*",err,error))// &
              & ") is not equal to the number of dynamic matrices ("// &
              & trim(numbertovstring(dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES,"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Equations matrices dynamic matrices is not associated.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF
    
    exits("EquationsMatrices_DynamicStructureTypeSet")
    RETURN
999 errorsexits("EquationsMatrices_DynamicStructureTypeSet",err,error)
    RETURN 1
    
  END SUBROUTINE equationsmatrices_dynamicstructuretypeset
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equationsmatrices_linearstructuretypeset(EQUATIONS_MATRICES,STRUCTURE_TYPE,ERR,ERROR,*)
    
    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(IN) :: STRUCTURE_TYPE(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
    TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    TYPE(varying_string) :: LOCAL_ERROR

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

    IF(ASSOCIATED(equations_matrices)) THEN
      IF(equations_matrices%EQUATIONS_MATRICES_FINISHED) THEN
        CALL flagerror("Equations matrices have been finished.",err,error,*999)
      ELSE
        linear_matrices=>equations_matrices%LINEAR_MATRICES
        IF(ASSOCIATED(linear_matrices)) THEN
          IF(SIZE(structure_type,1)==linear_matrices%NUMBER_OF_LINEAR_MATRICES) THEN
            DO matrix_idx=1,linear_matrices%NUMBER_OF_LINEAR_MATRICES
              equations_matrix=>linear_matrices%MATRICES(matrix_idx)%PTR
              IF(ASSOCIATED(equations_matrix)) THEN
                SELECT CASE(structure_type(matrix_idx))
                CASE(equations_matrix_no_structure)
                  equations_matrix%STRUCTURE_TYPE=equations_matrix_no_structure
                CASE(equations_matrix_fem_structure)
                  equations_matrix%STRUCTURE_TYPE=equations_matrix_fem_structure
                CASE(equations_matrix_diagonal_structure)
                  equations_matrix%STRUCTURE_TYPE=equations_matrix_diagonal_structure
                CASE(equations_matrix_nodal_structure)
                  equations_matrix%STRUCTURE_TYPE=equations_matrix_nodal_structure
                CASE DEFAULT
                  local_error="The specified strucutre type of "// &
                    & trim(numbertovstring(structure_type(matrix_idx),"*",err,error))//" for linear matrix number "// &
                    & trim(numbertovstring(matrix_idx,"*",err,error))//" is invalid."
                  CALL flagerror(local_error,err,error,*999)
                END SELECT
              ELSE
                CALL flagerror("Equations matrix is not associated.",err,error,*999)
              ENDIF
            ENDDO !matrix_idx
          ELSE
            local_error="The size of the structure type array ("//trim(numbertovstring(SIZE(structure_type,1),"*",err,error))// &
              & ") is not equal to the number of linear matrices ("// &
              & trim(numbertovstring(linear_matrices%NUMBER_OF_LINEAR_MATRICES,"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Equations matrices linear matrices is not associated.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF
    
    exits("EquationsMatrices_LinearStructureTypeSet")
    RETURN
999 errorsexits("EquationsMatrices_LinearStructureTypeSet",err,error)
    RETURN 1
    
  END SUBROUTINE equationsmatrices_linearstructuretypeset
  
  !
  !================================================================================================================================
  !

  SUBROUTINE equationsmatrices_nonlinearstructuretypeset0(EQUATIONS_MATRICES,STRUCTURE_TYPE,ERR,ERROR,*)
    
    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(IN) :: STRUCTURE_TYPE(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
    TYPE(equations_jacobian_type), POINTER :: JACOBIAN_MATRIX
    TYPE(varying_string) :: LOCAL_ERROR

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

    IF(ASSOCIATED(equations_matrices)) THEN
      IF(equations_matrices%EQUATIONS_MATRICES_FINISHED) THEN
        CALL flagerror("Equations matrices have been finished.",err,error,*999)
      ELSE
        nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
        IF(ASSOCIATED(nonlinear_matrices)) THEN
          IF(SIZE(structure_type,1)==nonlinear_matrices%NUMBER_OF_JACOBIANS) THEN
            DO matrix_idx=1,nonlinear_matrices%NUMBER_OF_JACOBIANS
              jacobian_matrix=>nonlinear_matrices%JACOBIANS(matrix_idx)%PTR
              IF(ASSOCIATED(jacobian_matrix)) THEN
                SELECT CASE(structure_type(matrix_idx))
                CASE(equations_matrix_no_structure)
                  jacobian_matrix%STRUCTURE_TYPE=equations_matrix_no_structure
                CASE(equations_matrix_fem_structure)
                  jacobian_matrix%STRUCTURE_TYPE=equations_matrix_fem_structure
                CASE(equations_matrix_diagonal_structure)
                  jacobian_matrix%STRUCTURE_TYPE=equations_matrix_diagonal_structure
                CASE(equations_matrix_nodal_structure)
                  jacobian_matrix%STRUCTURE_TYPE=equations_matrix_nodal_structure
                CASE DEFAULT
                  local_error="The specified strucutre type of "// &
                    & trim(numbertovstring(structure_type(matrix_idx),"*",err,error))//" for the Jacobian matrix is invalid."
                  CALL flagerror(local_error,err,error,*999)
                END SELECT
              ELSE
                CALL flagerror("Equations matrix is not associated.",err,error,*999)
              ENDIF
            ENDDO
          ELSE
            local_error="The size of the structure type array ("//trim(numbertovstring(SIZE(structure_type,1),"*",err,error))// &
              & ") is not equal to the number of Jacobian matrices ("// &
              & trim(numbertovstring(nonlinear_matrices%NUMBER_OF_JACOBIANS,"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Equations matrices nonlinear matrices is not associated.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF
    
    exits("EquationsMatrices_NonlinearStructureTypeSet0")
    RETURN
999 errors("EquationsMatrices_NonlinearStructureTypeSet0",err,error)
    exits("EquationsMatrices_NonlinearStructureTypeSet0")
    RETURN 1
    
  END SUBROUTINE equationsmatrices_nonlinearstructuretypeset0

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

  SUBROUTINE equationsmatrices_nonlinearstructuretypeset1(EQUATIONS_MATRICES,STRUCTURE_TYPE,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(IN) :: STRUCTURE_TYPE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG), ALLOCATABLE :: STRUCTURE_TYPES(:)
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES

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

    IF(ASSOCIATED(equations_matrices)) THEN
      IF(equations_matrices%EQUATIONS_MATRICES_FINISHED) THEN
        CALL flagerror("Equations matrices have been finished.",err,error,*999)
      ELSE
        nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
        IF(ASSOCIATED(nonlinear_matrices)) THEN
          ALLOCATE(structure_types(nonlinear_matrices%NUMBER_OF_JACOBIANS),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate storage types.",err,error,*999)
          structure_types=structure_type
          CALL equationsmatrices_nonlinearstructuretypeset0(equations_matrices,structure_types,err,error,*999)
          DEALLOCATE(structure_types)
        ELSE
          CALL flagerror("Equations matrices nonlinear matrices is not associated.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF

    exits("EquationsMatrices_NonlinearStructureTypeSet1")
    RETURN
999 errors("EquationsMatrices_NonlinearStructureTypeSet1",err,error)
    exits("EquationsMatrices_NonlinearStructureTypeSet1")
    RETURN 1
    
  END SUBROUTINE equationsmatrices_nonlinearstructuretypeset1

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

  SUBROUTINE equations_matrices_finalise(EQUATIONS_MATRICES,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
   
    enters("EQUATIONS_MATRICES_FINALISE",err,error,*999)

    IF(ASSOCIATED(equations_matrices)) THEN
      CALL equations_matrices_dynamic_finalise(equations_matrices%DYNAMIC_MATRICES,err,error,*999)
      CALL equations_matrices_linear_finalise(equations_matrices%LINEAR_MATRICES,err,error,*999)
      CALL equations_matrices_nonlinear_finalise(equations_matrices%NONLINEAR_MATRICES,err,error,*999)
      CALL equations_matrices_rhs_finalise(equations_matrices%RHS_VECTOR,err,error,*999)      
      CALL equations_matrices_source_finalise(equations_matrices%SOURCE_VECTOR,err,error,*999)      
      DEALLOCATE(equations_matrices)
    ENDIF
       
    exits("EQUATIONS_MATRICES_FINALISE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_finalise

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

  SUBROUTINE equations_matrices_initialise(EQUATIONS,ERR,ERROR,*)
    
     !Argument variables
    TYPE(equations_type), POINTER :: EQUATIONS
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DUMMY_ERR
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(varying_string) :: DUMMY_ERROR
    
    enters("EQUATIONS_MATRICES_INITIALISE",err,error,*998)

    IF(ASSOCIATED(equations)) THEN
      IF(ASSOCIATED(equations%EQUATIONS_MATRICES)) THEN
        CALL flagerror("Equations matrices is already associated for this equations.",err,error,*998)
      ELSE
        equations_mapping=>equations%EQUATIONS_MAPPING
        IF(ASSOCIATED(equations_mapping)) THEN
          IF(equations_mapping%EQUATIONS_MAPPING_FINISHED) THEN
            ALLOCATE(equations%EQUATIONS_MATRICES,stat=err)
            IF(err/=0) CALL flagerror("Could not allocate equations equations matrices.",err,error,*999)
            equations%EQUATIONS_MATRICES%EQUATIONS=>equations
            equations%EQUATIONS_MATRICES%EQUATIONS_MATRICES_FINISHED=.false.
            equations%EQUATIONS_MATRICES%EQUATIONS_MAPPING=>equations_mapping
            NULLIFY(equations%EQUATIONS_MATRICES%SOLVER_MAPPING)
            equations%EQUATIONS_MATRICES%NUMBER_OF_ROWS=equations_mapping%NUMBER_OF_ROWS
            equations%EQUATIONS_MATRICES%TOTAL_NUMBER_OF_ROWS=equations_mapping%TOTAL_NUMBER_OF_ROWS
            equations%EQUATIONS_MATRICES%NUMBER_OF_GLOBAL_ROWS=equations_mapping%NUMBER_OF_GLOBAL_ROWS
            NULLIFY(equations%EQUATIONS_MATRICES%DYNAMIC_MATRICES)
            NULLIFY(equations%EQUATIONS_MATRICES%LINEAR_MATRICES)
            NULLIFY(equations%EQUATIONS_MATRICES%NONLINEAR_MATRICES)
            NULLIFY(equations%EQUATIONS_MATRICES%RHS_VECTOR)
            NULLIFY(equations%EQUATIONS_MATRICES%SOURCE_VECTOR)            
            CALL equations_matrices_dynamic_initialise(equations%EQUATIONS_MATRICES,err,error,*999)            
            CALL equations_matrices_linear_initialise(equations%EQUATIONS_MATRICES,err,error,*999)            
            CALL equations_matrices_nonlinear_initialise(equations%EQUATIONS_MATRICES,err,error,*999)            
            CALL equations_matrices_rhs_initialise(equations%EQUATIONS_MATRICES,err,error,*999)            
            CALL equations_matrices_source_initialise(equations%EQUATIONS_MATRICES,err,error,*999)            
          ELSE
            CALL flagerror("Equations mapping has not been finished.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Equations equations mapping is not associated.",err,error,*998)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations is not associated.",err,error,*998)
    ENDIF
       
    exits("EQUATIONS_MATRICES_INITIALISE")
    RETURN
999 CALL equations_matrices_finalise(equations%EQUATIONS_MATRICES,dummy_err,dummy_error,*998)
998 errorsexits("EQUATIONS_MATRICES_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_initialise

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

  SUBROUTINE equations_matrices_values_initialise(EQUATIONS_MATRICES,SELECTION_TYPE,VALUE,ERR,ERROR,*)
    
    !Argument variables
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    INTEGER(INTG), INTENT(IN) :: SELECTION_TYPE
    REAL(DP), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: matrix_idx
    TYPE(equations_jacobian_type), POINTER :: JACOBIAN_MATRIX
    TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
    TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONLINEAR_MATRICES
    TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
    TYPE(equations_matrices_source_type), POINTER :: SOURCE_VECTOR
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    
    enters("EQUATIONS_MATRICES_VALUES_INITIALISE",err,error,*999)
    
    IF(ASSOCIATED(equations_matrices)) THEN
      IF(selection_type==equations_matrices_all.OR. &
        & selection_type==equations_matrices_dynamic_only.OR. &
        & selection_type==equations_matrices_linear_only.OR. &
        & selection_type==equations_matrices_nonlinear_only) THEN
        dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
        IF(ASSOCIATED(dynamic_matrices)) THEN
          DO matrix_idx=1,dynamic_matrices%NUMBER_OF_DYNAMIC_MATRICES
            equations_matrix=>dynamic_matrices%MATRICES(matrix_idx)%PTR
            IF(ASSOCIATED(equations_matrix)) THEN
              IF(equations_matrix%UPDATE_MATRIX) THEN
                CALL distributed_matrix_all_values_set(equations_matrix%MATRIX,VALUE,err,error,*999)
              ENDIF
            ELSE
              CALL flagerror("Equations matrix is not associated.",err,error,*999)
            ENDIF
          ENDDO !matrix_idx
        ENDIF
      ENDIF
      IF(selection_type==equations_matrices_all.OR. &
        & selection_type==equations_matrices_dynamic_only.OR. &
        & selection_type==equations_matrices_linear_only.OR. &
        & selection_type==equations_matrices_nonlinear_only) THEN
        linear_matrices=>equations_matrices%LINEAR_MATRICES
        IF(ASSOCIATED(linear_matrices)) THEN
          DO matrix_idx=1,linear_matrices%NUMBER_OF_LINEAR_MATRICES
            equations_matrix=>linear_matrices%MATRICES(matrix_idx)%PTR
            IF(ASSOCIATED(equations_matrix)) THEN
              IF(equations_matrix%UPDATE_MATRIX) THEN
                CALL distributed_matrix_all_values_set(equations_matrix%MATRIX,VALUE,err,error,*999)
              ENDIF
            ELSE
              CALL flagerror("Equations matrix is not associated.",err,error,*999)
            ENDIF
          ENDDO !matrix_idx
        ENDIF
      ENDIF
      IF(selection_type==equations_matrices_all.OR. &
        & selection_type==equations_matrices_nonlinear_only.OR. &
        & selection_type==equations_matrices_jacobian_only) THEN
        nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
        IF(ASSOCIATED(nonlinear_matrices)) THEN
          DO matrix_idx=1,nonlinear_matrices%NUMBER_OF_JACOBIANS
            jacobian_matrix=>nonlinear_matrices%JACOBIANS(matrix_idx)%PTR
            IF(ASSOCIATED(jacobian_matrix)) THEN
              IF(jacobian_matrix%UPDATE_JACOBIAN) THEN
                CALL distributed_matrix_all_values_set(jacobian_matrix%JACOBIAN,VALUE,err,error,*999)
              ENDIF
            ELSE
              CALL flagerror("Jacobian matrix is not associated.",err,error,*999)
            ENDIF
          ENDDO
        ENDIF
      ENDIF
      IF(selection_type==equations_matrices_all.OR. &
        & selection_type==equations_matrices_nonlinear_only.OR. &
        & selection_type==equations_matrices_residual_only.OR. &
        & selection_type==equations_matrices_rhs_residual_only.OR. &
        & selection_type==equations_matrices_residual_source_only.OR. &
        & selection_type==equations_matrices_vectors_only) THEN
        nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
        IF(ASSOCIATED(nonlinear_matrices)) THEN
          IF(nonlinear_matrices%UPDATE_RESIDUAL) THEN
            CALL distributed_vector_all_values_set(nonlinear_matrices%RESIDUAL,VALUE,err,error,*999)
          ENDIF
        ENDIF
      ENDIF
      IF(selection_type==equations_matrices_all.OR. &
        & selection_type==equations_matrices_dynamic_only.OR. &
        & selection_type==equations_matrices_linear_only.OR. &
        & selection_type==equations_matrices_nonlinear_only.OR. &
        & selection_type==equations_matrices_rhs_only.OR. &
        & selection_type==equations_matrices_rhs_residual_only.OR. &
        & selection_type==equations_matrices_rhs_source_only.OR. &
        & selection_type==equations_matrices_vectors_only) THEN    
        rhs_vector=>equations_matrices%RHS_VECTOR
        IF(ASSOCIATED(rhs_vector)) THEN
          IF(rhs_vector%UPDATE_VECTOR) THEN
            CALL distributed_vector_all_values_set(rhs_vector%VECTOR,VALUE,err,error,*999)
          ENDIF
        ENDIF
      ENDIF
      IF(selection_type==equations_matrices_all.OR. &
        & selection_type==equations_matrices_dynamic_only.OR. &
        & selection_type==equations_matrices_linear_only.OR. &
        & selection_type==equations_matrices_nonlinear_only.OR. &
        & selection_type==equations_matrices_source_only.OR. &
        & selection_type==equations_matrices_rhs_source_only.OR. &
        & selection_type==equations_matrices_residual_source_only.OR. &
        & selection_type==equations_matrices_vectors_only) THEN    
        source_vector=>equations_matrices%SOURCE_VECTOR
        IF(ASSOCIATED(source_vector)) THEN
          IF(source_vector%UPDATE_VECTOR) THEN
            CALL distributed_vector_all_values_set(source_vector%VECTOR,VALUE,err,error,*999)
          ENDIF
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Equations matrices is not associated.",err,error,*999)
    ENDIF

    exits("EQUATIONS_MATRICES_VALUES_INITIALISE")
    RETURN
999 errorsexits("EQUATIONS_MATRICES_VALUES_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE equations_matrices_values_initialise

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

  SUBROUTINE equationsmatrix_structurecalculate(equationsMatrix,numberOfNonZeros,rowIndices,columnIndices,list,err,error,*)

    !Argument variables
    TYPE(equations_matrix_type), POINTER :: equationsMatrix
    INTEGER(INTG), INTENT(OUT) :: numberOfNonZeros
    INTEGER(INTG), POINTER :: rowIndices(:)
    INTEGER(INTG), POINTER :: columnIndices(:)
    type(linkedlist),pointer :: list(:) 
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) ::  columnIdx,dummyErr,elementIdx,globalColumn,localColumn,local_ny,matrixNumber,mk,mp,ne,nh,nh2,nn,nnk,np
    INTEGER(INTG) ::  numberOfColumns,nyy,nyyg,npg,nhg,local_cols,local_dof,mv
    INTEGER(INTG) :: dofIdx,nodeIdx,componentIdx,localDofIdx
    INTEGER(INTG) :: versionIdx,derivativeIdx,numberOfDerivatives,numberOfVersions
    INTEGER(INTG), ALLOCATABLE :: columns(:)
    REAL(DP) :: sparsity
    TYPE(basis_type), POINTER :: basis
    TYPE(domain_mapping_type), POINTER :: dependentDofsDomainMapping
    TYPE(domain_elements_type), POINTER :: domainElements
    TYPE(domain_nodes_type), POINTER :: domainNodes
    TYPE(equations_type), POINTER :: equations
    TYPE(equations_mapping_type), POINTER :: equationsMapping
    TYPE(equations_mapping_dynamic_type), POINTER :: dynamicMapping
    TYPE(equations_mapping_linear_type), POINTER :: linearMapping
    TYPE(equations_matrices_type), POINTER :: equationsMatrices
    TYPE(equations_matrices_dynamic_type), POINTER :: dynamicMatrices
    TYPE(equations_matrices_linear_type), POINTER :: linearMatrices
    TYPE(equations_set_type), POINTER :: equationsSet
    TYPE(field_type), POINTER :: dependentField
    TYPE(field_dof_to_param_map_type), POINTER :: dependentDofsParamMapping
    TYPE(field_variable_type), POINTER :: fieldVariable
    TYPE(list_ptr_type), ALLOCATABLE :: columnIndicesLists(:)
    TYPE(varying_string) :: dummyError,localError
    enters("EquationsMatrix_StructureCalculate",err,error,*998)

    numberofnonzeros=0
    IF(ASSOCIATED(equationsmatrix)) THEN
      IF(.NOT.ASSOCIATED(rowindices)) THEN
        IF(.NOT.ASSOCIATED(columnindices)) THEN
          matrixnumber=equationsmatrix%MATRIX_NUMBER
          SELECT CASE(equationsmatrix%STRUCTURE_TYPE)
          CASE(equations_matrix_no_structure)
            CALL flagerror("There is no structure to calculate for a matrix with no structure.",err,error,*998)
          CASE(equations_matrix_fem_structure)
            SELECT CASE(equationsmatrix%STORAGE_TYPE)
            CASE(matrix_compressed_row_storage_type)
              linearmatrices=>equationsmatrix%LINEAR_MATRICES
              dynamicmatrices=>equationsmatrix%DYNAMIC_MATRICES
              IF(ASSOCIATED(dynamicmatrices).OR.ASSOCIATED(linearmatrices)) THEN
                IF(ASSOCIATED(dynamicmatrices)) THEN
                  equationsmatrices=>dynamicmatrices%EQUATIONS_MATRICES
                ELSE
                  equationsmatrices=>linearmatrices%EQUATIONS_MATRICES
                ENDIF
                IF(ASSOCIATED(equationsmatrices)) THEN
                  equations=>equationsmatrices%EQUATIONS
                  IF(ASSOCIATED(equations)) THEN
                    equationsmapping=>equationsmatrices%EQUATIONS_MAPPING
                    IF(ASSOCIATED(equationsmapping)) THEN
                      dynamicmapping=>equationsmapping%DYNAMIC_MAPPING
                      linearmapping=>equationsmapping%LINEAR_MAPPING
                      IF(ASSOCIATED(dynamicmapping).OR.ASSOCIATED(linearmapping)) THEN
                        equationsset=>equations%EQUATIONS_SET
                        IF(ASSOCIATED(equationsset)) THEN
                          dependentfield=>equationsset%DEPENDENT%DEPENDENT_FIELD
                          IF(ASSOCIATED(dependentfield)) THEN
                            IF(ASSOCIATED(dynamicmatrices)) THEN
                              fieldvariable=>dynamicmapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrixnumber)%VARIABLE
                            ELSE
                              fieldvariable=>linearmapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrixnumber)%VARIABLE
                            ENDIF
                            IF(ASSOCIATED(fieldvariable)) THEN
                              dependentdofsdomainmapping=>fieldvariable%DOMAIN_MAPPING
                              IF(ASSOCIATED(dependentdofsdomainmapping)) THEN
                                dependentdofsparammapping=>fieldvariable%DOF_TO_PARAM_MAP
                                IF(ASSOCIATED(dependentdofsparammapping)) THEN
                                  !Allocate lists
                                  ALLOCATE(columnindiceslists(dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL),stat=err)
                                  IF(err/=0) CALL flagerror("Could not allocate column indices lists.",err,error,*999)
                                  !Allocate row indices
                                  ALLOCATE(rowindices(dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL+1),stat=err)
                                  IF(err/=0) CALL flagerror("Could not allocate row indices.",err,error,*999)
                                  rowindices(1)=1
                                  
                                  !First, loop over the rows and calculate the number of non-zeros
                                  numberofnonzeros=0
                                  DO local_ny=1,dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL
                                    IF(dependentdofsparammapping%DOF_TYPE(1,local_ny)==field_node_dof_type) THEN
                                      nyy=dependentdofsparammapping%DOF_TYPE(2,local_ny)!value for a particular field dof (local_ny)
                                      np=dependentdofsparammapping%NODE_DOF2PARAM_MAP(3,nyy)!node number (np) of the field parameter
                                      nh=dependentdofsparammapping%NODE_DOF2PARAM_MAP(4,nyy)!component number (nh) of the field parameter
                                      domainnodes=>fieldvariable%COMPONENTS(nh)%DOMAIN%TOPOLOGY%NODES
                                      
                                      !Set up list
                                      NULLIFY(columnindiceslists(local_ny)%PTR)
                                      CALL list_create_start(columnindiceslists(local_ny)%PTR,err,error,*999)
                                      CALL list_data_type_set(columnindiceslists(local_ny)%PTR,list_intg_type,err,error,*999)
                                      CALL list_initial_size_set(columnindiceslists(local_ny)%PTR,domainnodes%NODES(np)% &
                                        & number_of_surrounding_elements*fieldvariable%COMPONENTS(nh)% &
                                        & maxnumberelementinterpolationparameters,err,error,*999)
                                      CALL list_create_finish(columnindiceslists(local_ny)%PTR,err,error,*999)
                                      !Loop over all elements containing the dof
                                      DO elementidx=1,domainnodes%NODES(np)%NUMBER_OF_SURROUNDING_ELEMENTS
                                        ne=domainnodes%NODES(np)%SURROUNDING_ELEMENTS(elementidx)
                                        DO nh2=1,fieldvariable%NUMBER_OF_COMPONENTS
                                          domainelements=>fieldvariable%COMPONENTS(nh2)%DOMAIN%TOPOLOGY%ELEMENTS
                                          basis=>domainelements%ELEMENTS(ne)%BASIS
                                          DO nn=1,basis%NUMBER_OF_NODES
                                            mp=domainelements%ELEMENTS(ne)%ELEMENT_NODES(nn)
                                            DO nnk=1,basis%NUMBER_OF_DERIVATIVES(nn)
                                              mk=domainelements%ELEMENTS(ne)%ELEMENT_DERIVATIVES(nnk,nn)
                                              mv=domainelements%ELEMENTS(ne)%elementVersions(nnk,nn)
                                              !Find the local and global column and add the global column to the indices list
                                              localcolumn=fieldvariable%COMPONENTS(nh2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% &
                                                & nodes(mp)%DERIVATIVES(mk)%VERSIONS(mv)
                                              globalcolumn=fieldvariable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(localcolumn)
                                          
                                              CALL list_item_add(columnindiceslists(local_ny)%PTR,globalcolumn,err,error,*999)
                                                
                                            ENDDO !mk
                                          ENDDO !nn
                                        ENDDO !nh2
                                      ENDDO !elementIdx
                                      CALL list_remove_duplicates(columnindiceslists(local_ny)%PTR,err,error,*999)
                                      CALL list_number_of_items_get(columnindiceslists(local_ny)%PTR,numberofcolumns, &
                                        & err,error,*999)
                                      numberofnonzeros=numberofnonzeros+numberofcolumns
                                      rowindices(local_ny+1)=numberofnonzeros+1
                                    ELSE
                                      localerror="Local dof number "//trim(numbertovstring(local_ny,"*",err,error))// &
                                        & " is not a node based dof."
                                      CALL flagerror(localerror,err,error,*999)
                                    ENDIF
                                  ENDDO !local_ny
                                  
                                  
                                  !Allocate and setup the column locations
                                  ALLOCATE(columnindices(numberofnonzeros),stat=err)

                                  ALLOCATE(list(dependentdofsdomainmapping%NUMBER_OF_GLOBAL))

                                  IF(err/=0) CALL flagerror("Could not allocate column indices.",err,error,*999)
                                  DO local_ny=1,dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL
                                   
                                    CALL list_detach_and_destroy(columnindiceslists(local_ny)%PTR,numberofcolumns,columns, &
                                      & err,error,*999)        
                                    DO columnidx=1,numberofcolumns
                                      !COLUMNS store the list of nonzero column indices for each local row (local_ny)
                                      columnindices(rowindices(local_ny)+columnidx-1)=columns(columnidx) 

                                      ! global to local columns
                                       IF(ASSOCIATED(linearmapping).OR.ASSOCIATED(dynamicmapping)) THEN
                                         IF(ASSOCIATED(dynamicmatrices)) THEN
                                           local_cols=equationsmatrices%equations_mapping%dynamic_mapping &
                                             & %equations_matrix_to_var_maps(1)%column_dofs_mapping%global_to_local_map &
                                             & (columns(columnidx))%LOCAL_NUMBER(1)
                                           local_dof = local_cols
                                           ! Column to dof mapping?
                                           !local_dof=equationsMatrices%equations_mapping%dynamic_mapping% &
                                            ! & equations_matrix_to_var_maps(1)%column_to_dof_map(local_cols)
                                         ELSE
                                           local_cols=equationsmatrices%equations_mapping%linear_mapping &
                                             & %equations_matrix_to_var_maps(1)%column_dofs_mapping%global_to_local_map &
                                             & (columns(columnidx))%LOCAL_NUMBER(1)
                                           local_dof = local_cols
                                         ENDIF
                                       ENDIF
                                       nyyg=dependentdofsparammapping%DOF_TYPE(2,local_dof)
                                       npg=dependentdofsparammapping%NODE_DOF2PARAM_MAP(3,nyyg)
                                       nhg=dependentdofsparammapping%NODE_DOF2PARAM_MAP(4,nyyg)
                                       domainnodes=>fieldvariable%COMPONENTS(nhg)%DOMAIN%TOPOLOGY%NODES
                            
                                      ! Check whether boundary node    
                                      IF(domainnodes%NODES(npg)%BOUNDARY_NODE)THEN
                                        CALL linkedlist_add(list(columns(columnidx)),local_ny,err,error,*999)
                                      ENDIF
                                    
                                    ENDDO !columnIdx
                                    DEALLOCATE(columns)                                    
                                  ENDDO !local_ny

                                 
                                  IF(diagnostics1) THEN
                                    CALL write_string(diagnostic_output_type,"Equations matrix structure:",err,error,*999)
                                    CALL write_string_value(diagnostic_output_type,"Equations matrix number : ",matrixnumber, &
                                      & err,error,*999)
                                    CALL write_string_value(diagnostic_output_type,"  Number of rows = ", &
                                      & dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL,err,error,*999)
                                    CALL write_string_value(diagnostic_output_type,"  Number of columns = ", &
                                      & dependentdofsdomainmapping%NUMBER_OF_GLOBAL,err,error,*999)
                                    CALL write_string_value(diagnostic_output_type,"  Number of non zeros = ", &
                                      & numberofnonzeros,err,error,*999)
                                    IF(dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL* &
                                      & dependentdofsdomainmapping%NUMBER_OF_GLOBAL/=0) THEN
                                      sparsity=(1.0_dp-REAL(numberofnonzeros,dp)/REAL(dependentdofsdomainmapping% &
                                        & TOTAL_NUMBER_OF_LOCAL*dependentDofsDomainMapping%NUMBER_OF_GLOBAL,DP))*100.0_DP
                                      CALL WRITE_STRING_FMT_VALUE(diagnostic_output_type,"  Sparsity (% of zeros) = ", &
                                        & sparsity,"F6.2",err,error,*999)
                                    ENDIF
                                    CALL write_string_vector(diagnostic_output_type,1,1,dependentdofsdomainmapping% &
                                      & total_number_of_local+1,8,8,rowindices,'("  Row indices    :",8(X,I13))', &
                                      & '(18X,8(X,I13))',err,error,*999)
                                    CALL write_string_vector(diagnostic_output_type,1,1,numberofnonzeros,8,8,columnindices, &
                                      & '("  Column indices :",8(X,I13))','(18X,8(X,I13))', err,error,*999)
                                  ENDIF
                                ELSE
                                  CALL flagerror("Dependent dofs parameter mapping is not associated.",err,error,*999)
                                ENDIF
                              ELSE
                                CALL flagerror("Dependent dofs domain mapping is not associated.",err,error,*999)
                              ENDIF
                            ELSE
                              CALL flagerror("Dependent field variable is not associated.",err,error,*999)
                            ENDIF
                          ELSE
                            CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
                          ENDIF
                        ELSE
                          CALL flagerror("Equations set is not associated.",err,error,*999)
                        ENDIF
                      ELSE
                        CALL flagerror("Either equations mapping dynamic mapping or linear mapping is not associated.", &
                          & err,error,*999)
                      ENDIF
                    ELSE
                      CALL flagerror("Equations mapping is not associated.",err,error,*999)
                    ENDIF
                  ELSE
                    CALL flagerror("Equations is not associated.",err,error,*999)
                  ENDIF
                ELSE
                  CALL flagerror("Dynamic or linear matrices equations matrices is not associated.",err,error,*999)
                ENDIF
              ELSE
                CALL flagerror("Either equations matrix dynamic or linear matrices is not associated.",err,error,*999)
              ENDIF
            CASE DEFAULT
              localerror="The matrix storage type of "// &
                & trim(numbertovstring(equationsmatrix%STORAGE_TYPE,"*",err,error))//" is invalid."
              CALL flagerror(localerror,err,error,*999)
            END SELECT

          CASE(equations_matrix_nodal_structure)
            SELECT CASE(equationsmatrix%STORAGE_TYPE)
            CASE(matrix_compressed_row_storage_type)
              linearmatrices=>equationsmatrix%LINEAR_MATRICES
              dynamicmatrices=>equationsmatrix%DYNAMIC_MATRICES
              IF(ASSOCIATED(dynamicmatrices).OR.ASSOCIATED(linearmatrices)) THEN
                IF(ASSOCIATED(dynamicmatrices)) THEN
                  equationsmatrices=>dynamicmatrices%EQUATIONS_MATRICES
                ELSE
                  equationsmatrices=>linearmatrices%EQUATIONS_MATRICES
                ENDIF
                IF(ASSOCIATED(equationsmatrices)) THEN
                  equations=>equationsmatrices%EQUATIONS
                  IF(ASSOCIATED(equations)) THEN
                    equationsmapping=>equationsmatrices%EQUATIONS_MAPPING
                    IF(ASSOCIATED(equationsmapping)) THEN
                      dynamicmapping=>equationsmapping%DYNAMIC_MAPPING
                      linearmapping=>equationsmapping%LINEAR_MAPPING
                      IF(ASSOCIATED(dynamicmapping).OR.ASSOCIATED(linearmapping)) THEN
                        equationsset=>equations%EQUATIONS_SET
                        IF(ASSOCIATED(equationsset)) THEN
                          dependentfield=>equationsset%DEPENDENT%DEPENDENT_FIELD
                          IF(ASSOCIATED(dependentfield)) THEN
                            IF(ASSOCIATED(dynamicmatrices)) THEN
                              fieldvariable=>dynamicmapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrixnumber)%VARIABLE
                            ELSE
                              fieldvariable=>linearmapping%EQUATIONS_MATRIX_TO_VAR_MAPS(matrixnumber)%VARIABLE
                            ENDIF
                            IF(ASSOCIATED(fieldvariable)) THEN
                              dependentdofsdomainmapping=>fieldvariable%DOMAIN_MAPPING
                              IF(ASSOCIATED(dependentdofsdomainmapping)) THEN
                                dependentdofsparammapping=>fieldvariable%DOF_TO_PARAM_MAP
                                IF(ASSOCIATED(dependentdofsparammapping)) THEN
                                  !Allocate lists
                                  ALLOCATE(columnindiceslists(dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL),stat=err)
                                  IF(err/=0) CALL flagerror("Could not allocate column indices lists.",err,error,*999)
                                  !Allocate row indices
                                  ALLOCATE(rowindices(dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL+1),stat=err)
                                  IF(err/=0) CALL flagerror("Could not allocate row indices.",err,error,*999)
                                  rowindices(1)=1
                                  
                                  !First, loop over the rows and calculate the number of non-zeros
                                  numberofnonzeros=0
                                  DO localdofidx=1,dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL
                                    IF(dependentdofsparammapping%DOF_TYPE(1,localdofidx)==field_node_dof_type) THEN
                                      dofidx=dependentdofsparammapping%DOF_TYPE(2,localdofidx)!value for a particular field dof (localDofIdx)
                                      nodeidx=dependentdofsparammapping%NODE_DOF2PARAM_MAP(3,dofidx)!node number (np) of the field parameter
                                      componentidx=dependentdofsparammapping%NODE_DOF2PARAM_MAP(4,dofidx)!component number (nh) of the field parameter
                                      domainnodes=>fieldvariable%COMPONENTS(componentidx)%DOMAIN%TOPOLOGY%NODES
                                      
                                      !Set up list
                                      NULLIFY(columnindiceslists(localdofidx)%PTR)
                                      CALL list_create_start(columnindiceslists(localdofidx)%PTR,err,error,*999)
                                      CALL list_data_type_set(columnindiceslists(localdofidx)%PTR,list_intg_type,err,error,*999)

                                      CALL list_initial_size_set(columnindiceslists(localdofidx)%PTR, &
                                        & fieldvariable%NUMBER_OF_COMPONENTS* &
                                        & fieldvariable%maxNumberElementInterpolationParameters,err,error,*999)

                                      CALL list_create_finish(columnindiceslists(localdofidx)%PTR,err,error,*999)
                                      !Loop over all components,nodes,derivatives, and versions
                                      DO componentidx=1,fieldvariable%NUMBER_OF_COMPONENTS
                                        numberofderivatives=fieldvariable%components(componentidx)%domain%topology% &
                                         & nodes%nodes(nodeidx)%NUMBER_OF_DERIVATIVES
                                        DO derivativeidx=1,numberofderivatives
                                          numberofversions=fieldvariable%components(componentidx)%domain%topology% &
                                           & nodes%nodes(nodeidx)%derivatives(derivativeidx)%numberOfVersions
                                          DO versionidx=1,numberofversions
                                            localcolumn=fieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP% &
                                             & node_param2dof_map%NODES(nodeidx)%DERIVATIVES(derivativeidx)%VERSIONS(versionidx)
                                            globalcolumn=fieldvariable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(localcolumn)

                                            CALL list_item_add(columnindiceslists(localdofidx)%PTR,globalcolumn,err,error,*999)
                                            
                                          ENDDO !versionIdx
                                        ENDDO !derivativeIdx
                                      ENDDO !componentIdx

                                      CALL list_remove_duplicates(columnindiceslists(localdofidx)%PTR,err,error,*999)
                                      CALL list_number_of_items_get(columnindiceslists(localdofidx)%PTR,numberofcolumns, &
                                        & err,error,*999)
                                      numberofnonzeros=numberofnonzeros+numberofcolumns
                                      rowindices(localdofidx+1)=numberofnonzeros+1
                                    ELSE
                                      localerror="Local dof number "//trim(numbertovstring(localdofidx,"*",err,error))// &
                                        & " is not a node based dof."
                                      CALL flagerror(localerror,err,error,*999)
                                    ENDIF
                                  ENDDO !localDofIdx
                                  
                                  !Allocate and setup the column locations
                                  ALLOCATE(columnindices(numberofnonzeros),stat=err)
                                  ALLOCATE(list(dependentdofsdomainmapping%NUMBER_OF_GLOBAL))

                                  IF(err/=0) CALL flagerror("Could not allocate column indices.",err,error,*999)
                                  DO localdofidx=1,dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL
                                   
                                    CALL list_detach_and_destroy(columnindiceslists(localdofidx)%PTR,numberofcolumns,columns, &
                                      & err,error,*999)        
                                    DO columnidx=1,numberofcolumns
                                      !columns store the list of nonzero column indices for each local row (localDofIdx)
                                      columnindices(rowindices(localdofidx)+columnidx-1)=columns(columnidx) 

                                      ! global to local columns
                                       IF(ASSOCIATED(linearmapping).OR.ASSOCIATED(dynamicmapping)) THEN
                                         IF(ASSOCIATED(dynamicmatrices)) THEN
                                           local_cols=equationsmatrices%equations_mapping%dynamic_mapping &
                                             & %equations_matrix_to_var_maps(1)%column_dofs_mapping%global_to_local_map &
                                             & (columns(columnidx))%LOCAL_NUMBER(1)
                                           local_dof = local_cols
                                           ! Column to dof mapping?
                                           !local_dof=equationsMatrices%equations_mapping%dynamic_mapping% &
                                            ! & equations_matrix_to_var_maps(1)%column_to_dof_map(local_cols)
                                         ELSE
                                           local_cols=equationsmatrices%equations_mapping%linear_mapping &
                                             & %equations_matrix_to_var_maps(1)%column_dofs_mapping%global_to_local_map &
                                             & (columns(columnidx))%LOCAL_NUMBER(1)
                                           local_dof = local_cols
                                         ENDIF
                                       ENDIF
                                       nyyg=dependentdofsparammapping%DOF_TYPE(2,local_dof)
                                       npg=dependentdofsparammapping%NODE_DOF2PARAM_MAP(3,nyyg)
                                       nhg=dependentdofsparammapping%NODE_DOF2PARAM_MAP(4,nyyg)
                                       domainnodes=>fieldvariable%COMPONENTS(nhg)%DOMAIN%TOPOLOGY%NODES
                            
                                      ! Check whether boundary node    
                                      IF(domainnodes%NODES(npg)%BOUNDARY_NODE)THEN
                                        CALL linkedlist_add(list(columns(columnidx)),localdofidx,err,error,*999)
                                      ENDIF
                                    
                                    ENDDO !columnIdx
                                    DEALLOCATE(columns)                                    
                                  ENDDO !localDofIdx
                                 
                                  IF(diagnostics1) THEN
                                    CALL write_string(diagnostic_output_type,"Equations matrix structure:",err,error,*999)
                                    CALL write_string_value(diagnostic_output_type,"Equations matrix number : ",matrixnumber, &
                                      & err,error,*999)
                                    CALL write_string_value(diagnostic_output_type,"  Number of rows = ", &
                                      & dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL,err,error,*999)
                                    CALL write_string_value(diagnostic_output_type,"  Number of columns = ", &
                                      & dependentdofsdomainmapping%NUMBER_OF_GLOBAL,err,error,*999)
                                    CALL write_string_value(diagnostic_output_type,"  Number of non zeros = ", &
                                      & numberofnonzeros,err,error,*999)
                                    IF(dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL* &
                                      & dependentdofsdomainmapping%NUMBER_OF_GLOBAL/=0) THEN
                                      sparsity=(1.0_dp-REAL(numberofnonzeros,dp)/REAL(dependentdofsdomainmapping% &
                                        & TOTAL_NUMBER_OF_LOCAL*dependentDofsDomainMapping%NUMBER_OF_GLOBAL,DP))*100.0_DP
                                      CALL WRITE_STRING_FMT_VALUE(diagnostic_output_type,"  Sparsity (% of zeros) = ", &
                                        & sparsity,"F6.2",err,error,*999)
                                    ENDIF
                                    CALL write_string_vector(diagnostic_output_type,1,1,dependentdofsdomainmapping% &
                                      & total_number_of_local+1,8,8,rowindices,'("  Row indices    :",8(X,I13))', &
                                      & '(18X,8(X,I13))',err,error,*999)
                                    CALL write_string_vector(diagnostic_output_type,1,1,numberofnonzeros,8,8,columnindices, &
                                      & '("  Column indices :",8(X,I13))','(18X,8(X,I13))', err,error,*999)
                                  ENDIF
                                ELSE
                                  CALL flagerror("Dependent dofs parameter mapping is not associated.",err,error,*999)
                                ENDIF
                              ELSE
                                CALL flagerror("Dependent dofs domain mapping is not associated.",err,error,*999)
                              ENDIF
                            ELSE
                              CALL flagerror("Dependent field variable is not associated.",err,error,*999)
                            ENDIF
                          ELSE
                            CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
                          ENDIF
                        ELSE
                          CALL flagerror("Equations set is not associated.",err,error,*999)
                        ENDIF
                      ELSE
                        CALL flagerror("Either equations mapping dynamic mapping or linear mapping is not associated.", &
                          & err,error,*999)
                      ENDIF
                    ELSE
                      CALL flagerror("Equations mapping is not associated.",err,error,*999)
                    ENDIF
                  ELSE
                    CALL flagerror("Equations is not associated.",err,error,*999)
                  ENDIF
                ELSE
                  CALL flagerror("Dynamic or linear matrices equations matrices is not associated.",err,error,*999)
                ENDIF
              ELSE
                CALL flagerror("Either equations matrix dynamic or linear matrices is not associated.",err,error,*999)
              ENDIF

            CASE DEFAULT
              localerror="The matrix storage type of "// &
                & trim(numbertovstring(equationsmatrix%STORAGE_TYPE,"*",err,error))//" is invalid."
              CALL flagerror(localerror,err,error,*999)
            END SELECT

          CASE(equations_matrix_diagonal_structure)
            CALL flagerror("There is not structure to calculate for a diagonal matrix.",err,error,*998)
          CASE DEFAULT
            localerror="The matrix structure type of "// &
              & trim(numbertovstring(equationsmatrix%STRUCTURE_TYPE,"*",err,error))//" is invalid."
            CALL flagerror(localerror,err,error,*998)
          END SELECT
        ELSE
          CALL flagerror("Column indices is already associated.",err,error,*998)
        ENDIF
      ELSE
        CALL flagerror("Row indieces is already associated.",err,error,*998)
      ENDIF
    ELSE
      CALL flagerror("Equations matrix is not associated.",err,error,*999)
    ENDIF
      
    exits("EquationsMatrix_StructureCalculate")
    RETURN
999 IF(ASSOCIATED(rowindices)) DEALLOCATE(rowindices)
    IF(ASSOCIATED(columnindices)) DEALLOCATE(columnindices)
    IF(ALLOCATED(columns)) DEALLOCATE(columns)
    IF(ALLOCATED(columnindiceslists)) THEN
      DO localdofidx=1,dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL
        IF(ASSOCIATED(columnindiceslists(localdofidx)%PTR)) &
          & CALL list_destroy(columnindiceslists(localdofidx)%PTR,dummyerr,dummyerror,*998)
      ENDDO !localDofIdx
      DEALLOCATE(columnindiceslists)
    ENDIF
998 errorsexits("EquationsMatrix_StructureCalculate",err,error)
    RETURN 1
  END SUBROUTINE equationsmatrix_structurecalculate

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

  SUBROUTINE jacobianmatrix_structurecalculate(jacobianMatrix,numberOfNonZeros,rowIndices,columnIndices,err,error,*)

    !Argument variables
    TYPE(equations_jacobian_type), POINTER :: jacobianMatrix
    INTEGER(INTG), INTENT(OUT) :: numberOfNonZeros
    INTEGER(INTG), POINTER :: rowIndices(:)
    INTEGER(INTG), POINTER :: columnIndices(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) ::  columnIdx,dummyErr,elementIdx,globalColumn,localColumn,local_ny,mk,mp,ne,nh,nh2,nn,nnk,np,mv, &
      & numberOfColumns,nyy,matrixNumber
    INTEGER(INTG) :: dofIdx,nodeIdx,componentIdx,versionIdx,derivativeIdx,numberOfVersions,numberOfDerivatives
    INTEGER(INTG) :: localDofIdx
    INTEGER(INTG), ALLOCATABLE :: columns(:)
    REAL(DP) :: sparsity
    TYPE(basis_type), POINTER :: basis
    TYPE(domain_mapping_type), POINTER :: dependentDofsDomainMapping,rowDofsDomainMapping
    TYPE(domain_elements_type), POINTER :: domainElements
    TYPE(domain_nodes_type), POINTER :: domainNodes
    TYPE(equations_type), POINTER :: equations
    TYPE(equations_mapping_type), POINTER :: equationsMapping
    TYPE(equations_mapping_nonlinear_type), POINTER :: nonlinearMapping
    TYPE(equations_matrices_type), POINTER :: equationsMatrices
    TYPE(equations_matrices_nonlinear_type), POINTER :: NONlinearMatrices
    TYPE(equations_set_type), POINTER :: equationsSet
    TYPE(field_type), POINTER :: dependentField
    TYPE(field_dof_to_param_map_type), POINTER :: dependentDofsParamMapping,rowDofsParamMapping
    TYPE(field_variable_type), POINTER :: fieldVariable,rowVariable
    TYPE(list_ptr_type), ALLOCATABLE :: columnIndicesLists(:)
    TYPE(varying_string) :: dummyError,localError

    enters("JacobianMatrix_StructureCalculate",err,error,*998)

    numberofnonzeros=0
    IF(ASSOCIATED(jacobianmatrix)) THEN
      matrixnumber=jacobianmatrix%JACOBIAN_NUMBER
      IF(.NOT.ASSOCIATED(rowindices)) THEN
        IF(.NOT.ASSOCIATED(columnindices)) THEN
          SELECT CASE(jacobianmatrix%STRUCTURE_TYPE)
          CASE(equations_matrix_no_structure)
            CALL flagerror("Not implemented.",err,error,*998)
          CASE(equations_matrix_fem_structure)
            SELECT CASE(jacobianmatrix%STORAGE_TYPE)
            CASE(matrix_compressed_row_storage_type)
              nonlinearmatrices=>jacobianmatrix%NONLINEAR_MATRICES
              IF(ASSOCIATED(nonlinearmatrices)) THEN
                equationsmatrices=>nonlinearmatrices%EQUATIONS_MATRICES
                IF(ASSOCIATED(equationsmatrices)) THEN
                  equations=>equationsmatrices%EQUATIONS
                  IF(ASSOCIATED(equations)) THEN
                    equationsmapping=>equationsmatrices%EQUATIONS_MAPPING
                    IF(ASSOCIATED(equationsmapping)) THEN
                      nonlinearmapping=>equationsmapping%NONLINEAR_MAPPING
                      IF(ASSOCIATED(nonlinearmapping)) THEN
                        equationsset=>equations%EQUATIONS_SET
                        IF(ASSOCIATED(equationsset)) THEN
                          dependentfield=>equationsset%DEPENDENT%DEPENDENT_FIELD
                          IF(ASSOCIATED(dependentfield)) THEN
                            fieldvariable=>nonlinearmapping%JACOBIAN_TO_VAR_MAP(matrixnumber)%VARIABLE
                            IF(ASSOCIATED(fieldvariable)) THEN
                              dependentdofsdomainmapping=>fieldvariable%DOMAIN_MAPPING
                              IF(ASSOCIATED(dependentdofsdomainmapping)) THEN
                                dependentdofsparammapping=>fieldvariable%DOF_TO_PARAM_MAP
                                IF(ASSOCIATED(dependentdofsparammapping)) THEN
                                  !If RHS variable exists, use this for row DOFs, else use the first nonlinear variable
                                  IF(ASSOCIATED(equationsmapping%RHS_MAPPING)) THEN
                                    rowvariable=>equationsmapping%RHS_MAPPING%RHS_VARIABLE
                                  ELSE
                                    rowvariable=>nonlinearmapping%JACOBIAN_TO_VAR_MAP(1)%VARIABLE
                                  ENDIF
                                  IF(ASSOCIATED(rowvariable)) THEN
                                    rowdofsdomainmapping=>rowvariable%DOMAIN_MAPPING
                                    rowdofsparammapping=>rowvariable%DOF_TO_PARAM_MAP
                                  ELSE
                                    CALL flagerror("RHS or first nonlinear variable is not associated",err,error,*999)
                                  ENDIF
                                  IF(ASSOCIATED(rowdofsdomainmapping)) THEN
                                    IF(ASSOCIATED(rowdofsparammapping)) THEN
                                      !Allocate lists
                                      ALLOCATE(columnindiceslists(rowdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL),stat=err)
                                      IF(err/=0) CALL flagerror("Could not allocate column indices lists.",err,error,*999)
                                      !Allocate row indices
                                      ALLOCATE(rowindices(rowdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL+1),stat=err)
                                      IF(err/=0) CALL flagerror("Could not allocate row indices.",err,error,*999)
                                      rowindices(1)=1
                                      !First, loop over the rows and calculate the number of non-zeros
                                      numberofnonzeros=0
                                      DO local_ny=1,rowdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL
                                        SELECT CASE(rowdofsparammapping%DOF_TYPE(1,local_ny))
                                        CASE(field_constant_interpolation)
                                          CALL flagerror("Constant interpolation is not implemented yet.",err,error,*999)
                                        CASE(field_node_dof_type)
                                          nyy=rowdofsparammapping%DOF_TYPE(2,local_ny)
                                          np=rowdofsparammapping%NODE_DOF2PARAM_MAP(3,nyy) !node number
                                          nh=rowdofsparammapping%NODE_DOF2PARAM_MAP(4,nyy) !component number
                                          domainnodes=>rowvariable%COMPONENTS(nh)%DOMAIN%TOPOLOGY%NODES
                                          !Set up list
                                          NULLIFY(columnindiceslists(local_ny)%PTR)
                                          CALL list_create_start(columnindiceslists(local_ny)%PTR,err,error,*999)
                                          CALL list_data_type_set(columnindiceslists(local_ny)%PTR,list_intg_type,err,error,*999)
                                          CALL list_initial_size_set(columnindiceslists(local_ny)%PTR,domainnodes%NODES(np)% &
                                            & number_of_surrounding_elements*rowvariable%COMPONENTS(nh)% &
                                            & maxnumberelementinterpolationparameters,err,error,*999)
                                          CALL list_create_finish(columnindiceslists(local_ny)%PTR,err,error,*999)
                                          !Loop over all elements containing the dof
                                          DO elementidx=1,domainnodes%NODES(np)%NUMBER_OF_SURROUNDING_ELEMENTS
                                            ne=domainnodes%NODES(np)%SURROUNDING_ELEMENTS(elementidx)
                                            DO nh2=1,fieldvariable%NUMBER_OF_COMPONENTS
                                              SELECT CASE(fieldvariable%COMPONENTS(nh2)%INTERPOLATION_TYPE)
                                              CASE(field_constant_interpolation)
                                                ! do nothing? this will probably never be encountered...?
                                              CASE(field_element_based_interpolation)
                                                localcolumn=fieldvariable%COMPONENTS(nh2)%PARAM_TO_DOF_MAP% &
                                                  & element_param2dof_map%ELEMENTS(ne)
                                                globalcolumn=fieldvariable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(localcolumn)
                                                CALL list_item_add(columnindiceslists(local_ny)%PTR,globalcolumn,err,error,*999)
                                              CASE(field_node_based_interpolation)
                                                domainelements=>fieldvariable%COMPONENTS(nh2)%DOMAIN%TOPOLOGY%ELEMENTS
                                                basis=>domainelements%ELEMENTS(ne)%BASIS
                                                DO nn=1,basis%NUMBER_OF_NODES
                                                  mp=domainelements%ELEMENTS(ne)%ELEMENT_NODES(nn)
                                                  DO nnk=1,basis%NUMBER_OF_DERIVATIVES(nn)
                                                    mk=domainelements%ELEMENTS(ne)%ELEMENT_DERIVATIVES(nnk,nn)
                                                    mv=domainelements%ELEMENTS(ne)%elementVersions(nnk,nn)
                                                    !Find the local and global column and add the global column to the indices list
                                                    localcolumn=fieldvariable%COMPONENTS(nh2)%PARAM_TO_DOF_MAP% &
                                                      & node_param2dof_map%NODES(mp)%DERIVATIVES(mk)%VERSIONS(mv)
                                                    globalcolumn=fieldvariable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(localcolumn)
                                                    CALL list_item_add(columnindiceslists(local_ny)%PTR,globalcolumn, &
                                                      & err,error,*999)
                                                  ENDDO !mk
                                                ENDDO !nn
                                              CASE(field_grid_point_based_interpolation)
                                                CALL flagerror("Grid point based interpolation is not implemented yet.",& 
                                                  & err,error,*999)
                                              CASE(field_gauss_point_based_interpolation)
                                                CALL flagerror("Gauss point based interpolation is not implemented yet.",&
                                                  & err,error,*999)
                                              CASE DEFAULT
                                                localerror="Local dof number "//trim(numbertovstring(local_ny,"*",err,error))// &
                                                  & " has invalid interpolation type."
                                                CALL flagerror(localerror,err,error,*999)
                                              END SELECT
                                            ENDDO !nh2
                                          ENDDO !elementIdx
                                          CALL list_remove_duplicates(columnindiceslists(local_ny)%PTR,err,error,*999)
                                          CALL list_number_of_items_get(columnindiceslists(local_ny)%PTR,numberofcolumns, &
                                            & err,error,*999)
                                          numberofnonzeros=numberofnonzeros+numberofcolumns
                                          rowindices(local_ny+1)=numberofnonzeros+1
                                        CASE(field_element_dof_type)
                                          ! row corresponds to a variable that's element-wisely interpolated
                                          nyy=rowdofsparammapping%DOF_TYPE(2,local_ny)          ! nyy = index in ELEMENT_DOF2PARAM_MAP
                                          ne=rowdofsparammapping%ELEMENT_DOF2PARAM_MAP(1,nyy)   ! current element (i.e. corresponds to current dof)
                                          nh=rowdofsparammapping%ELEMENT_DOF2PARAM_MAP(2,nyy)   ! current variable component
                                          domainelements=>rowvariable%COMPONENTS(nh)%DOMAIN%TOPOLOGY%ELEMENTS
                                          basis=>domainelements%ELEMENTS(ne)%BASIS
                                          !Set up list
                                          NULLIFY(columnindiceslists(local_ny)%PTR)
                                          CALL list_create_start(columnindiceslists(local_ny)%PTR,err,error,*999)
                                          CALL list_data_type_set(columnindiceslists(local_ny)%PTR,list_intg_type,err,error,*999)
                                          CALL list_initial_size_set(columnindiceslists(local_ny)%PTR, &
                                            & rowvariable%COMPONENTS(nh)%maxNumberElementInterpolationParameters+1, &
                                            & err,error,*999) ! size = all nodal dofs + itself
                                          CALL list_create_finish(columnindiceslists(local_ny)%PTR,err,error,*999)
                                          DO nh2=1,fieldvariable%NUMBER_OF_COMPONENTS
                                            SELECT CASE(fieldvariable%COMPONENTS(nh2)%INTERPOLATION_TYPE)
                                            CASE(field_constant_interpolation)
                                              CALL flagerror("Constant interpolation is not implemented yet.",err,error,*999)
                                            CASE(field_element_based_interpolation)
                                              ! it's assumed that element-based variables arne't directly coupled
                                              ! put a diagonal entry
                                              localcolumn=fieldvariable%COMPONENTS(nh2)%PARAM_TO_DOF_MAP% &
                                                & element_param2dof_map%ELEMENTS(ne)
                                              globalcolumn=fieldvariable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(localcolumn)
                                              CALL list_item_add(columnindiceslists(local_ny)%PTR,globalcolumn,err,error,*999)
                                            CASE(field_node_based_interpolation)
                                              ! loop over all nodes in the element (and dofs belonging to them)
                                              DO nn=1,basis%NUMBER_OF_NODES
                                                mp=domainelements%ELEMENTS(ne)%ELEMENT_NODES(nn)
                                                DO nnk=1,basis%NUMBER_OF_DERIVATIVES(nn)
                                                  mk=domainelements%ELEMENTS(ne)%ELEMENT_DERIVATIVES(nnk,nn)
                                                  mv=domainelements%ELEMENTS(ne)%elementVersions(nnk,nn)
                                                  !Find the local and global column and add the global column to the indices list
                                                  localcolumn=fieldvariable%COMPONENTS(nh2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% &
                                                    & nodes(mp)%DERIVATIVES(mk)%VERSIONS(mv)
                                                  globalcolumn=fieldvariable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(localcolumn)
                                                  CALL list_item_add(columnindiceslists(local_ny)%PTR,globalcolumn, &
                                                    & err,error,*999)
                                                ENDDO !mk
                                              ENDDO !nn
                                            CASE(field_grid_point_based_interpolation)
                                              CALL flagerror("Grid point based interpolation is not implemented yet.", &
                                                & err,error,*999)
                                            CASE(field_gauss_point_based_interpolation)
                                              CALL flagerror("Gauss point based interpolation is not implemented yet.", &
                                                & err,error,*999)
                                            CASE DEFAULT
                                              localerror="Local dof number "//trim(numbertovstring(local_ny,"*",err,error))// &
                                                & " has invalid interpolation type."
                                              CALL flagerror(localerror,err,error,*999)
                                            END SELECT
                                          ENDDO !nh2
                                          ! clean up the list
                                          CALL list_remove_duplicates(columnindiceslists(local_ny)%PTR,err,error,*999)
                                          CALL list_number_of_items_get(columnindiceslists(local_ny)%PTR,numberofcolumns, &
                                            & err,error,*999)
                                          numberofnonzeros=numberofnonzeros+numberofcolumns
                                          rowindices(local_ny+1)=numberofnonzeros+1
                                        CASE(field_grid_point_based_interpolation)
                                          CALL flagerror("Grid point based interpolation is not implemented yet.",err,error,*999)
                                        CASE(field_gauss_point_based_interpolation)
                                          CALL flagerror("Gauss point based interpolation is not implemented yet.",err,error,*999)
                                        CASE DEFAULT
                                          localerror="Local dof number "//trim(numbertovstring(local_ny,"*",err,error))// &
                                            & " has an invalid type."
                                          CALL flagerror(localerror,err,error,*999)
                                        END SELECT
                                      ENDDO !local_ny
                                      !Allocate and setup the column locations
                                      ALLOCATE(columnindices(numberofnonzeros),stat=err)
                                      IF(err/=0) CALL flagerror("Could not allocate column indices.",err,error,*999)
                                      DO local_ny=1,rowdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL
                                        CALL list_detach_and_destroy(columnindiceslists(local_ny)%PTR,numberofcolumns,columns, &
                                          & err,error,*999)
                                        DO columnidx=1,numberofcolumns
                                          columnindices(rowindices(local_ny)+columnidx-1)=columns(columnidx)
                                        ENDDO !columnIdx
                                        DEALLOCATE(columns)
                                      ENDDO !local_ny
                                      IF(diagnostics1) THEN
                                        CALL write_string(diagnostic_output_type,"Jacobian matrix structure:",err,error,*999)
                                        CALL write_string_value(diagnostic_output_type,"  Number of rows = ", &
                                          & dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL,err,error,*999)
                                        CALL write_string_value(diagnostic_output_type,"  Number of columns = ", &
                                          & dependentdofsdomainmapping%NUMBER_OF_GLOBAL,err,error,*999)
                                        CALL write_string_value(diagnostic_output_type,"  Number of non zeros = ", &
                                          & numberofnonzeros,err,error,*999)
                                        IF(dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL* &
                                          & dependentdofsdomainmapping%NUMBER_OF_GLOBAL/=0) THEN
                                          sparsity=(1.0_dp-REAL(numberofnonzeros,dp)/REAL(dependentdofsdomainmapping% &
                                            & TOTAL_NUMBER_OF_LOCAL*dependentDofsDomainMapping%NUMBER_OF_GLOBAL,DP))*100.0_DP
                                          CALL WRITE_STRING_FMT_VALUE(diagnostic_output_type,"  Sparsity (% of zeros) = ", &
                                            & sparsity,"F6.2",err,error,*999)
                                        ENDIF
                                        CALL write_string_vector(diagnostic_output_type,1,1,rowdofsdomainmapping% &
                                          & total_number_of_local+1,8,8,rowindices,'("  Row indices    :",8(X,I13))', &
                                          & '(18X,8(X,I13))',err,error,*999)
                                        CALL write_string_vector(diagnostic_output_type,1,1,numberofnonzeros,8,8,columnindices,&
                                          & '("  Column indices :",8(X,I13))','(18X,8(X,I13))', err,error,*999)
                                      ENDIF
                                    ELSE
                                      CALL flagerror("Row dofs parameter mapping is not associated.",err,error,*999)
                                    ENDIF
                                  ELSE
                                    CALL flagerror("Row dofs domain mapping is not associated.",err,error,*999)
                                  ENDIF
                                ELSE
                                  CALL flagerror("Dependent dofs parameter mapping is not associated.",err,error,*999)
                                ENDIF
                              ELSE
                                CALL flagerror("Dependent dofs domain mapping is not associated.",err,error,*999)
                              ENDIF
                            ELSE
                              CALL flagerror("Dependent field variable is not associated.",err,error,*999)
                            ENDIF
                          ELSE
                            CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
                          ENDIF
                        ELSE
                          CALL flagerror("Equations set is not associated.",err,error,*999)
                        ENDIF
                      ELSE
                        CALL flagerror("Equations mapping nonlinear mapping is not associated.",err,error,*999)
                      ENDIF
                    ELSE
                      CALL flagerror("Equations mapping is not associated.",err,error,*999)
                    ENDIF
                  ELSE
                    CALL flagerror("Equations is not associated.",err,error,*999)
                  ENDIF
                ELSE
                  CALL flagerror("Nonlinear matrices equations matrices is not associated.",err,error,*999)
                ENDIF
              ELSE
                CALL flagerror("Equations matrix nonlinear matrices is not associated.",err,error,*999)
              ENDIF
            CASE DEFAULT
              localerror="The matrix storage type of "// &
                & trim(numbertovstring(jacobianmatrix%STORAGE_TYPE,"*",err,error))//" is invalid."
              CALL flagerror(localerror,err,error,*999)
            END SELECT

          CASE(equations_matrix_nodal_structure)
            SELECT CASE(jacobianmatrix%STORAGE_TYPE)
            CASE(matrix_compressed_row_storage_type)
              nonlinearmatrices=>jacobianmatrix%NONLINEAR_MATRICES
              IF(ASSOCIATED(nonlinearmatrices)) THEN
                equationsmatrices=>nonlinearmatrices%EQUATIONS_MATRICES
                IF(ASSOCIATED(equationsmatrices)) THEN
                  equations=>equationsmatrices%EQUATIONS
                  IF(ASSOCIATED(equations)) THEN
                    equationsmapping=>equationsmatrices%EQUATIONS_MAPPING
                    IF(ASSOCIATED(equationsmapping)) THEN
                      nonlinearmapping=>equationsmapping%NONLINEAR_MAPPING
                      IF(ASSOCIATED(nonlinearmapping)) THEN
                        equationsset=>equations%EQUATIONS_SET
                        IF(ASSOCIATED(equationsset)) THEN
                          dependentfield=>equationsset%DEPENDENT%DEPENDENT_FIELD
                          IF(ASSOCIATED(dependentfield)) THEN
                            fieldvariable=>nonlinearmapping%JACOBIAN_TO_VAR_MAP(matrixnumber)%VARIABLE
                            IF(ASSOCIATED(fieldvariable)) THEN
                              dependentdofsdomainmapping=>fieldvariable%DOMAIN_MAPPING
                              IF(ASSOCIATED(dependentdofsdomainmapping)) THEN
                                dependentdofsparammapping=>fieldvariable%DOF_TO_PARAM_MAP
                                IF(ASSOCIATED(dependentdofsparammapping)) THEN
                                  !If RHS variable exists, use this for row DOFs, else use the first nonlinear variable
                                  IF(ASSOCIATED(equationsmapping%RHS_MAPPING)) THEN
                                    rowvariable=>equationsmapping%RHS_MAPPING%RHS_VARIABLE
                                  ELSE
                                    rowvariable=>nonlinearmapping%JACOBIAN_TO_VAR_MAP(1)%VARIABLE
                                  ENDIF
                                  IF(ASSOCIATED(rowvariable)) THEN
                                    rowdofsdomainmapping=>rowvariable%DOMAIN_MAPPING
                                    rowdofsparammapping=>rowvariable%DOF_TO_PARAM_MAP
                                  ELSE
                                    CALL flagerror("RHS or first nonlinear variable is not associated",err,error,*999)
                                  ENDIF
                                  IF(ASSOCIATED(rowdofsdomainmapping)) THEN
                                    IF(ASSOCIATED(rowdofsparammapping)) THEN
                                      !Allocate lists
                                      ALLOCATE(columnindiceslists(rowdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL),stat=err)
                                      IF(err/=0) CALL flagerror("Could not allocate column indices lists.",err,error,*999)
                                      !Allocate row indices
                                      ALLOCATE(rowindices(rowdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL+1),stat=err)
                                      IF(err/=0) CALL flagerror("Could not allocate row indices.",err,error,*999)
                                      rowindices(1)=1
                                      !First, loop over the rows and calculate the number of non-zeros
                                      numberofnonzeros=0
                                      DO localdofidx=1,rowdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL
                                        SELECT CASE(rowdofsparammapping%DOF_TYPE(1,localdofidx))
                                        CASE(field_constant_interpolation)
                                          CALL flagerror("Constant interpolation is not implemented yet.",err,error,*999)
                                        CASE(field_node_dof_type)
                                          dofidx=dependentdofsparammapping%DOF_TYPE(2,localdofidx)!value for a particular field dof (localDofIdx)
                                          nodeidx=dependentdofsparammapping%NODE_DOF2PARAM_MAP(3,dofidx)!node number (np) of the field parameter
                                          componentidx=dependentdofsparammapping%NODE_DOF2PARAM_MAP(4,dofidx)!component number (nh) of the field parameter
                                          domainnodes=>fieldvariable%COMPONENTS(componentidx)%DOMAIN%TOPOLOGY%NODES

                                          !Set up list
                                          NULLIFY(columnindiceslists(localdofidx)%PTR)
                                          CALL list_create_start(columnindiceslists(localdofidx)%PTR,err,error,*999)
                                          CALL list_data_type_set(columnindiceslists(localdofidx)%PTR,list_intg_type,err,error,*999)

                                          CALL list_initial_size_set(columnindiceslists(localdofidx)%PTR, &
                                            & fieldvariable%NUMBER_OF_COMPONENTS* &
                                            & fieldvariable%maxNumberElementInterpolationParameters,err,error,*999)

                                          CALL list_create_finish(columnindiceslists(localdofidx)%PTR,err,error,*999)
                                          !Loop over all components,nodes,derivatives, and versions
                                          DO componentidx=1,fieldvariable%NUMBER_OF_COMPONENTS
                                            SELECT CASE(fieldvariable%COMPONENTS(componentidx)%INTERPOLATION_TYPE)
                                            CASE(field_node_based_interpolation)
                                              numberofderivatives=fieldvariable%components(componentidx)%domain%topology% &
                                               & nodes%nodes(nodeidx)%NUMBER_OF_DERIVATIVES
                                              DO derivativeidx=1,numberofderivatives
                                                numberofversions=fieldvariable%components(componentidx)%domain%topology% &
                                                 & nodes%nodes(nodeidx)%derivatives(derivativeidx)%numberOfVersions
                                                DO versionidx=1,numberofversions
                                                  localcolumn=fieldvariable%COMPONENTS(componentidx)%PARAM_TO_DOF_MAP% &
                                                   & node_param2dof_map%NODES(nodeidx)%DERIVATIVES(derivativeidx)% &
                                                   & versions(versionidx)
                                                  globalcolumn=fieldvariable%DOMAIN_MAPPING%LOCAL_TO_GLOBAL_MAP(localcolumn)

                                                  CALL list_item_add(columnindiceslists(localdofidx)%PTR,globalcolumn, &
                                                    & err,error,*999)

                                                ENDDO !versionIdx
                                              ENDDO !derivativeIdx
                                            CASE DEFAULT
                                              localerror="Local dof number "//trim(numbertovstring(localdofidx,"*",err,error))// &
                                                & " has invalid interpolation type."
                                              CALL flagerror(localerror,err,error,*999)
                                            END SELECT
                                          ENDDO !componentIdx

                                          CALL list_remove_duplicates(columnindiceslists(localdofidx)%PTR,err,error,*999)
                                          CALL list_number_of_items_get(columnindiceslists(localdofidx)%PTR,numberofcolumns, &
                                            & err,error,*999)
                                          numberofnonzeros=numberofnonzeros+numberofcolumns
                                          rowindices(localdofidx+1)=numberofnonzeros+1
                                        CASE(field_element_dof_type)
                                          CALL flagerror("Element based interpolation is not implemented yet.",err,error,*999)
                                        CASE(field_grid_point_based_interpolation)
                                          CALL flagerror("Grid point based interpolation is not implemented yet.",err,error,*999)
                                        CASE(field_gauss_point_based_interpolation)
                                          CALL flagerror("Gauss point based interpolation is not implemented yet.",err,error,*999)
                                        CASE DEFAULT
                                          localerror="Local dof number "//trim(numbertovstring(localdofidx,"*",err,error))// &
                                            & " has an invalid type."
                                          CALL flagerror(localerror,err,error,*999)
                                        END SELECT
                                      ENDDO !localDofIdx
                                      !Allocate and setup the column locations
                                      ALLOCATE(columnindices(numberofnonzeros),stat=err)
                                      IF(err/=0) CALL flagerror("Could not allocate column indices.",err,error,*999)
                                      DO localdofidx=1,rowdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL
                                        CALL list_detach_and_destroy(columnindiceslists(localdofidx)%PTR,numberofcolumns,columns, &
                                          & err,error,*999)
                                        DO columnidx=1,numberofcolumns
                                          columnindices(rowindices(localdofidx)+columnidx-1)=columns(columnidx)
                                        ENDDO !columnIdx
                                        DEALLOCATE(columns)
                                      ENDDO !localDofIdx
                                      IF(diagnostics1) THEN
                                        CALL write_string(diagnostic_output_type,"Jacobian matrix structure :",err,error,*999)
                                        CALL write_string_value(diagnostic_output_type,"  Number of rows = ", &
                                          & dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL,err,error,*999)
                                        CALL write_string_value(diagnostic_output_type,"  Number of columns = ", &
                                          & dependentdofsdomainmapping%NUMBER_OF_GLOBAL,err,error,*999)
                                        CALL write_string_value(diagnostic_output_type,"  Number of non zeros = ", &
                                          & numberofnonzeros,err,error,*999)
                                        IF(dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL* &
                                          & dependentdofsdomainmapping%NUMBER_OF_GLOBAL/=0) THEN
                                          sparsity=(1.0_dp-REAL(numberofnonzeros,dp)/REAL(dependentdofsdomainmapping% &
                                            & TOTAL_NUMBER_OF_LOCAL*dependentDofsDomainMapping%NUMBER_OF_GLOBAL,DP))*100.0_DP
                                          CALL WRITE_STRING_FMT_VALUE(diagnostic_output_type,"  Sparsity (% of zeros) = ", &
                                            & sparsity,"F6.2",err,error,*999)
                                        ENDIF
                                        CALL write_string_vector(diagnostic_output_type,1,1,rowdofsdomainmapping% &
                                          & total_number_of_local+1,8,8,rowindices,'("  Row indices    :",8(X,I13))', &
                                          & '(18X,8(X,I13))',err,error,*999)
                                        CALL write_string_vector(diagnostic_output_type,1,1,numberofnonzeros,8,8,columnindices,&
                                          & '("  Column indices :",8(X,I13))','(18X,8(X,I13))', err,error,*999)
                                      ENDIF
                                    ELSE
                                      CALL flagerror("Row dofs parameter mapping is not associated.",err,error,*999)
                                    ENDIF
                                  ELSE
                                    CALL flagerror("Row dofs domain mapping is not associated.",err,error,*999)
                                  ENDIF
                                ELSE
                                  CALL flagerror("Dependent dofs parameter mapping is not associated.",err,error,*999)
                                ENDIF
                              ELSE
                                CALL flagerror("Dependent dofs domain mapping is not associated.",err,error,*999)
                              ENDIF
                            ELSE
                              CALL flagerror("Dependent field variable is not associated.",err,error,*999)
                            ENDIF
                          ELSE
                            CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
                          ENDIF
                        ELSE
                          CALL flagerror("Equations set is not associated.",err,error,*999)
                        ENDIF
                      ELSE
                        CALL flagerror("Equations mapping nonlinear mapping is not associated.",err,error,*999)
                      ENDIF
                    ELSE
                      CALL flagerror("Equations mapping is not associated.",err,error,*999)
                    ENDIF
                  ELSE
                    CALL flagerror("Equations is not associated.",err,error,*999)
                  ENDIF
                ELSE
                  CALL flagerror("Nonlinear matrices equations matrices is not associated.",err,error,*999)
                ENDIF
              ELSE
                CALL flagerror("Equations matrix nonlinear matrices is not associated.",err,error,*999)
              ENDIF
            CASE DEFAULT
              localerror="The matrix storage type of "// &
                & trim(numbertovstring(jacobianmatrix%STORAGE_TYPE,"*",err,error))//" is invalid."
              CALL flagerror(localerror,err,error,*999)
            END SELECT
          CASE DEFAULT
            localerror="The matrix structure type of "// &
              & trim(numbertovstring(jacobianmatrix%STRUCTURE_TYPE,"*",err,error))//" is invalid."
            CALL flagerror(localerror,err,error,*998)
          END SELECT
        ELSE
          CALL flagerror("Column indices is already associated.",err,error,*998)
        ENDIF
      ELSE
        CALL flagerror("Row indices is already associated.",err,error,*998)
      ENDIF
    ELSE
      CALL flagerror("Jacobian matrix is not associated.",err,error,*999)
    ENDIF
      
    exits("JacobianMatrix_StructureCalculate")
    RETURN
999 IF(ASSOCIATED(rowindices)) DEALLOCATE(rowindices)
    IF(ASSOCIATED(columnindices)) DEALLOCATE(columnindices)
    IF(ALLOCATED(columns)) DEALLOCATE(columns)
    IF(ALLOCATED(columnindiceslists)) THEN
      DO localdofidx=1,dependentdofsdomainmapping%TOTAL_NUMBER_OF_LOCAL
        IF(ASSOCIATED(columnindiceslists(localdofidx)%PTR)) &
          & CALL list_destroy(columnindiceslists(localdofidx)%PTR,dummyerr,dummyerror,*998)
      ENDDO !localDofIdx
      DEALLOCATE(columnindiceslists)
    ENDIF
998 errorsexits("JacobianMatrix_StructureCalculate",err,error)
    RETURN 1
  END SUBROUTINE jacobianmatrix_structurecalculate

  !
  !================================================================================================================================
  !
 
END MODULE equations_matrices_routines