analytic_analysis_routines.f90

Go to the documentation of this file.

MODULE analytic_analysis_routines

  USE basis_routines
  USE cmiss_mpi
  USE comp_environment
  USE constants
  USE field_routines
  USE input_output
  USE iso_varying_string
  USE kinds
  USE matrix_vector
#ifndef NOMPIMOD
  USE mpi
#endif
  USE strings
  USE timer
  USE types

#include "macros.h"

  IMPLICIT NONE

  PRIVATE

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

  !Module parameters

  INTEGER(INTG), PARAMETER :: absolute_error_type=1
  INTEGER(INTG), PARAMETER :: percentage_error_type=2
  INTEGER(INTG), PARAMETER :: relative_error_type=3

  !Module types

  !Module variables

  !Interfaces

  PUBLIC analyticanalysis_output
  
  PUBLIC analyticanalysis_absoluteerrorgetnode,analyticanalysis_percentageerrorgetnode, &
    & analyticanalysis_relativeerrorgetnode,analyticanalysis_rmserrorgetnode

  PUBLIC analyticanalysis_absoluteerrorgetelement,analyticanalysis_percentageerrorgetelement, &
    & analyticanalysis_relativeerrorgetelement,analyticanalysis_rmserrorgetelement

  PUBLIC analyticanalysis_absoluteerrorgetconstant,analyticanalysis_percentageerrorgetconstant, &
    & analyticanalysis_relativeerrorgetconstant
   
  PUBLIC analyticanalysis_integralnumericalvalueget,analyticanalysis_integralanalyticvalueget, &
    & analyticanalysis_integralpercentageerrorget,analyticanalysis_integralabsoluteerrorget, &
    & analyticanalysis_integralrelativeerrorget,analyticanalysis_integralnidnumericalvalueget, &
    & analyticanalysis_integralniderrorget
    
CONTAINS  

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

  SUBROUTINE analyticanalysis_output(FIELD,FILENAME,ERR,ERROR,*)
  
    !Argument variables 
    TYPE(field_type), INTENT(IN), POINTER :: FIELD
    CHARACTER(LEN=*) :: FILENAME
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: component_idx,deriv_idx,element_idx,GHOST_NUMBER(8),local_ny,MESH_COMPONENT,MPI_IERROR,node_idx, &
      & NUMBER(8),OUTPUT_ID,var_idx,variable_type
    REAL(DP) :: GHOST_RMS_ERROR_PER(8),GHOST_RMS_ERROR_ABS(8),GHOST_RMS_ERROR_REL(8),RMS_ERROR_PER(8),RMS_ERROR_ABS(8), &
      & RMS_ERROR_REL(8),VALUES(5)
    REAL(DP), POINTER :: ANALYTIC_VALUES(:),NUMERICAL_VALUES(:)
    REAL(DP), ALLOCATABLE :: INTEGRAL_ERRORS(:,:),GHOST_INTEGRAL_ERRORS(:,:)
    CHARACTER(LEN=40) :: FIRST_FORMAT
    CHARACTER(LEN=MAXSTRLEN) :: FILE_NAME
    TYPE(decomposition_type), POINTER :: DECOMPOSITION
    TYPE(decomposition_elements_type), POINTER :: ELEMENTS_DECOMPOSITION
    TYPE(decomposition_topology_type), POINTER :: DECOMPOSITION_TOPOLOGY
    TYPE(domain_type), POINTER :: DOMAIN
    TYPE(domain_elements_type), POINTER :: ELEMENTS_DOMAIN
    TYPE(domain_nodes_type), POINTER :: NODES_DOMAIN
    TYPE(domain_topology_type), POINTER :: DOMAIN_TOPOLOGY
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
    TYPE(varying_string) :: LOCAL_ERROR,LOCAL_STRING

    NULLIFY(analytic_values)
    NULLIFY(numerical_values)
    
    enters("AnalyticAnalysis_Output",err,error,*999)
    
    IF(ASSOCIATED(field)) THEN
      IF(field%FIELD_FINISHED) THEN
        IF(field%DEPENDENT_TYPE==field_dependent_type) THEN
          IF(len_trim(filename)>=1) THEN
!!TODO \todo have more general ascii file mechanism
            IF(computational_environment%NUMBER_COMPUTATIONAL_NODES>1) THEN
              WRITE(file_name,'(A,".opanal.",I0)') filename(1:len_trim(filename)),computational_environment% &
                & my_computational_node_number
            ELSE
              file_name=filename(1:len_trim(filename))//".opanal"
            ENDIF
            output_id=io1_file_unit
            OPEN(unit=output_id,file=file_name(1:len_trim(file_name)),status="REPLACE",form="FORMATTED",iostat=err)
            IF(err/=0) CALL flagerror("Error opening analysis output file.",err,error,*999)            
          ELSE
            output_id=general_output_type
          ENDIF
          decomposition=>field%DECOMPOSITION
          IF(ASSOCIATED(decomposition)) THEN
            decomposition_topology=>decomposition%TOPOLOGY
            IF(ASSOCIATED(decomposition_topology)) THEN
              CALL write_string(output_id,"Analytic error analysis:",err,error,*999)
              CALL write_string(output_id,"",err,error,*999)
              local_string="Field "//trim(number_to_vstring(field%USER_NUMBER,"*",err,error))//" : "//field%LABEL
              IF(err/=0) GOTO 999
              CALL write_string(output_id,local_string,err,error,*999)
              NULLIFY(numerical_values)
              NULLIFY(analytic_values)
              !Loop over the variables
              DO var_idx=1,field%NUMBER_OF_VARIABLES
                variable_type=field%VARIABLES(var_idx)%VARIABLE_TYPE
                field_variable=>field%VARIABLE_TYPE_MAP(variable_type)%PTR
                IF(ASSOCIATED(field_variable)) THEN
                  CALL write_string(output_id,"",err,error,*999)
                  local_string="Variable "//trim(number_to_vstring(variable_type,"*",err,error))//" : "// &
                    & field_variable%VARIABLE_LABEL
                  IF(err/=0) GOTO 999
                  CALL write_string(output_id,local_string,err,error,*999)
                  CALL write_string(output_id,"",err,error,*999)
                  !Get the dependent and analytic parameter sets
                  CALL field_parameter_set_data_get(field,variable_type,field_values_set_type,numerical_values,err,error,*999)
                  CALL field_parameter_set_data_get(field,variable_type,field_analytic_values_set_type,analytic_values, &
                    & err,error,*999)
                  !Loop over the components
                  DO component_idx=1,field%VARIABLES(var_idx)%NUMBER_OF_COMPONENTS
                    mesh_component=field_variable%COMPONENTS(component_idx)%MESH_COMPONENT_NUMBER
                    domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
                    IF(ASSOCIATED(domain)) THEN
                      domain_topology=>domain%TOPOLOGY
                      IF(ASSOCIATED(domain_topology)) THEN
                        local_string="Component "//trim(number_to_vstring(component_idx,"*",err,error))//" : "// &
                          & field_variable%COMPONENTS(component_idx)%COMPONENT_LABEL
                        IF(err/=0) GOTO 999
                        CALL write_string(output_id,local_string,err,error,*999)
                        CALL write_string(output_id,"",err,error,*999)
                        SELECT CASE(field%VARIABLES(var_idx)%COMPONENTS(component_idx)%INTERPOLATION_TYPE)
                        CASE(field_constant_interpolation)
                          CALL write_string(output_id,"Constant errors:",err,error,*999)
                          local_string="                       Numerical      Analytic       % error  Absolute err  Relative err"
                          CALL write_string(output_id,local_string,err,error,*999)
                          local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
                          values(1)=numerical_values(local_ny)
                          values(2)=analytic_values(local_ny)
                          values(3)=analytic_analysis_percentage_error(values(1),values(2))
                          values(4)=analytic_analysis_absolute_error(values(1),values(2))
                          values(5)=analytic_analysis_relative_error(values(1),values(2))
                          CALL write_string_vector(output_id,1,1,5,5,5,values,"(20X,5(2X,E12.5))","(20X,5(2X,E12.5))", &
                            & err,error,*999)
                        CASE(field_element_based_interpolation)
                          elements_domain=>domain_topology%ELEMENTS
                          IF(ASSOCIATED(elements_domain)) THEN
                            decomposition=>domain%DECOMPOSITION
                            IF(ASSOCIATED(decomposition)) THEN
                              decomposition_topology=>decomposition%TOPOLOGY
                              IF(ASSOCIATED(decomposition_topology)) THEN
                                elements_decomposition=>decomposition_topology%ELEMENTS
                                IF(ASSOCIATED(elements_decomposition)) THEN
                                  number=0
                                  rms_error_per=0.0_dp
                                  rms_error_abs=0.0_dp
                                  rms_error_rel=0.0_dp
                                  ghost_number=0
                                  ghost_rms_error_per=0.0_dp
                                  ghost_rms_error_abs=0.0_dp
                                  ghost_rms_error_rel=0.0_dp
                                  CALL write_string(output_id,"Element errors:",err,error,*999)
                                  local_string= &
                                    & "  Element#             Numerical      Analytic       % error  Absolute err  Relative err"
                                  CALL write_string(output_id,local_string,err,error,*999)
                                  DO element_idx=1,elements_domain%NUMBER_OF_ELEMENTS
                                    local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
                                      & element_param2dof_map%ELEMENTS(element_idx)
                                    values(1)=numerical_values(local_ny)
                                    values(2)=analytic_values(local_ny)
                                    values(3)=analytic_analysis_percentage_error(values(1),values(2))
                                    values(4)=analytic_analysis_absolute_error(values(1),values(2))
                                    values(5)=analytic_analysis_relative_error(values(1),values(2))
                                    !Accumlate the RMS errors
                                    number(1)=number(1)+1
                                    rms_error_per(1)=rms_error_per(1)+values(3)*values(3)
                                    rms_error_abs(1)=rms_error_abs(1)+values(4)*values(4)
                                    rms_error_rel(1)=rms_error_rel(1)+values(5)*values(5)
                                    WRITE(first_format,"(A,I10,A)") "('",elements_decomposition%ELEMENTS(element_idx)%USER_NUMBER, &
                                      & "',20X,3(2X,E12.5))"
                                    CALL write_string_vector(output_id,1,1,5,5,5,values,first_format,"(20X,5(2X,E12.5))", &
                                      & err,error,*999)
                                  ENDDO !element_idx
                                  DO element_idx=elements_domain%NUMBER_OF_ELEMENTS+1,elements_domain%TOTAL_NUMBER_OF_ELEMENTS
                                    local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
                                      & element_param2dof_map%ELEMENTS(element_idx)
                                    values(1)=numerical_values(local_ny)
                                    values(2)=analytic_values(local_ny)
                                    values(3)=analytic_analysis_percentage_error(values(1),values(2))
                                    values(4)=analytic_analysis_absolute_error(values(1),values(2))
                                    values(5)=analytic_analysis_relative_error(values(1),values(2))
                                    !Accumlate the RMS errors
                                    ghost_number(1)=ghost_number(1)+1
                                    ghost_rms_error_per(1)=ghost_rms_error_per(1)+values(3)*values(3)
                                    ghost_rms_error_abs(1)=ghost_rms_error_abs(1)+values(4)*values(4)
                                    ghost_rms_error_rel(1)=ghost_rms_error_rel(1)+values(5)*values(5)
                                    WRITE(first_format,"(A,I10,A)") "('",elements_decomposition%ELEMENTS(element_idx)%USER_NUMBER, &
                                      & "',20X,3(2X,E12.5))"
                                    CALL write_string_vector(output_id,1,1,5,5,5,values,first_format,"(20X,5(2X,E12.5))", &
                                      & err,error,*999)
                                  ENDDO !node_idx
                                  !Output RMS errors                  
                                  CALL write_string(output_id,"",err,error,*999)
                                  IF(number(1)>0) THEN
                                    IF(computational_environment%NUMBER_COMPUTATIONAL_NODES>1) THEN
                                      !Local elements only
                                      CALL write_string(output_id,"Local RMS errors:",err,error,*999)
                                      local_string= &
                                        & "                                                     % error  Absolute err  Relative err"
                                      CALL write_string(output_id,local_string,err,error,*999)
                                      values(1)=sqrt(rms_error_per(deriv_idx)/number(deriv_idx))
                                      values(2)=sqrt(rms_error_abs(deriv_idx)/number(deriv_idx))
                                      values(3)=sqrt(rms_error_rel(deriv_idx)/number(deriv_idx))
                                      CALL write_string_vector(output_id,1,1,3,3,3,values,"(46X,3(2X,E12.5))","(46X,3(2X,E12.5))", &
                                        & err,error,*999)
                                      !Local and ghost nodes
                                      CALL write_string(output_id,"Local + Ghost RMS errors:",err,error,*999)
                                      local_string= &
                                        & "                                                     % error  Absolute err  Relative err"
                                      CALL write_string(output_id,local_string,err,error,*999)
                                      values(1)=sqrt((rms_error_per(1)+ghost_rms_error_per(1))/(number(1)+ghost_number(1)))
                                      values(2)=sqrt((rms_error_abs(1)+ghost_rms_error_abs(1))/(number(1)+ghost_number(1)))
                                      values(3)=sqrt((rms_error_rel(1)+ghost_rms_error_rel(1))/(number(1)+ghost_number(1)))
                                      CALL write_string_vector(output_id,1,1,3,3,3,values,"(46X,3(2X,E12.5))","(46X,3(2X,E12.5))", &
                                        & err,error,*999)
                                      !Global RMS values
                                      !Collect the values across the ranks
                                      CALL mpi_allreduce(mpi_in_place,number,1,mpi_integer,mpi_sum, &
                                        & computational_environment%MPI_COMM,mpi_ierror)
                                      CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
                                      CALL mpi_allreduce(mpi_in_place,rms_error_per,1,mpi_double_precision,mpi_sum, &
                                        & computational_environment%MPI_COMM,mpi_ierror)
                                      CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
                                      CALL mpi_allreduce(mpi_in_place,rms_error_abs,1,mpi_double_precision,mpi_sum, &
                                        & computational_environment%MPI_COMM,mpi_ierror)
                                      CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
                                      CALL mpi_allreduce(mpi_in_place,rms_error_rel,1,mpi_double_precision,mpi_sum, &
                                        & computational_environment%MPI_COMM,mpi_ierror)
                                      CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
                                      CALL write_string(output_id,"Global RMS errors:",err,error,*999)
                                      local_string= &
                                        & "                                                     % error  Absolute err  Relative err"
                                      CALL write_string(output_id,local_string,err,error,*999)
                                      values(1)=sqrt(rms_error_per(1)/number(1))
                                      values(2)=sqrt(rms_error_abs(1)/number(1))
                                      values(3)=sqrt(rms_error_rel(1)/number(1))
                                      CALL write_string_vector(output_id,1,1,3,3,3,values,"(46X,3(2X,E12.5))","(46X,3(2X,E12.5))", &
                                        & err,error,*999)
                                    ELSE
                                      CALL write_string(output_id,"RMS errors:",err,error,*999)
                                      local_string= &
                                        & "                                                     % error  Absolute err  Relative err"
                                      CALL write_string(output_id,local_string,err,error,*999)
                                      values(1)=sqrt(rms_error_per(deriv_idx)/number(deriv_idx))
                                      values(2)=sqrt(rms_error_abs(deriv_idx)/number(deriv_idx))
                                      values(3)=sqrt(rms_error_rel(deriv_idx)/number(deriv_idx))
                                      CALL write_string_vector(output_id,1,1,3,3,3,values,"(46X,3(2X,E12.5))","(46X,3(2X,E12.5))", &
                                        & err,error,*999)
                                    ENDIF
                                  ENDIF
                                ELSE
                                  CALL flagerror("Decomposition topology elements is not associated.",err,error,*999)
                                ENDIF
                              ELSE
                                CALL flagerror("Decomposition topology is not associated.",err,error,*999)
                              ENDIF
                            ELSE
                              CALL flagerror("Domain decomposition is not associated.",err,error,*999)
                            ENDIF
                          ELSE
                            CALL flagerror("Elements domain topology is not associated.",err,error,*999)
                          ENDIF
                        CASE(field_node_based_interpolation)
                          nodes_domain=>domain_topology%NODES
                          IF(ASSOCIATED(nodes_domain)) THEN
                            number=0
                            rms_error_per=0.0_dp
                            rms_error_abs=0.0_dp
                            rms_error_rel=0.0_dp
                            ghost_number=0
                            ghost_rms_error_per=0.0_dp
                            ghost_rms_error_abs=0.0_dp
                            ghost_rms_error_rel=0.0_dp
                            CALL write_string(output_id,"Nodal errors:",err,error,*999)
                            local_string="     Node#  Deriv#     Numerical      Analytic       % error  Absolute err  Relative err"
                            CALL write_string(output_id,local_string,err,error,*999)
                            DO node_idx=1,nodes_domain%NUMBER_OF_NODES
                              DO deriv_idx=1,nodes_domain%NODES(node_idx)%NUMBER_OF_DERIVATIVES
                                local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
                                  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
                                values(1)=numerical_values(local_ny)
                                values(2)=analytic_values(local_ny)
                                values(3)=analytic_analysis_percentage_error(values(1),values(2))
                                values(4)=analytic_analysis_absolute_error(values(1),values(2))
                                values(5)=analytic_analysis_relative_error(values(1),values(2))
                                !Accumlate the RMS errors
                                number(deriv_idx)=number(deriv_idx)+1
                                rms_error_per(deriv_idx)=rms_error_per(deriv_idx)+values(3)*values(3)
                                rms_error_abs(deriv_idx)=rms_error_abs(deriv_idx)+values(4)*values(4)
                                rms_error_rel(deriv_idx)=rms_error_rel(deriv_idx)+values(5)*values(5)
                                IF(deriv_idx==1) THEN
                                  WRITE(first_format,"(A,I10,A,I6,A)") "('",nodes_domain%NODES(node_idx)%USER_NUMBER,"',2X,'", &
                                    & deriv_idx,"',5(2X,E12.5))"
                                ELSE
                                  WRITE(first_format,"(A,I6,A)") "(12X,'",deriv_idx,"',5(2X,E12.5))"
                                ENDIF
                                CALL write_string_vector(output_id,1,1,5,5,5,values,first_format,"(20X,5(2X,E12.5))",err,error,*999)
                              ENDDO !deriv_idx
                            ENDDO !node_idx
                            DO node_idx=nodes_domain%NUMBER_OF_NODES+1,nodes_domain%TOTAL_NUMBER_OF_NODES
                              DO deriv_idx=1,nodes_domain%NODES(node_idx)%NUMBER_OF_DERIVATIVES
                                local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
                                  & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
                                values(1)=numerical_values(local_ny)
                                values(2)=analytic_values(local_ny)
                                values(3)=analytic_analysis_percentage_error(values(1),values(2))
                                values(4)=analytic_analysis_absolute_error(values(1),values(2))
                                values(5)=analytic_analysis_relative_error(values(1),values(2))
                                !Accumlate the RMS errors
                                ghost_number(deriv_idx)=ghost_number(deriv_idx)+1
                                ghost_rms_error_per(deriv_idx)=ghost_rms_error_per(deriv_idx)+values(3)*values(3)
                                ghost_rms_error_abs(deriv_idx)=ghost_rms_error_abs(deriv_idx)+values(4)*values(4)
                                ghost_rms_error_rel(deriv_idx)=ghost_rms_error_rel(deriv_idx)+values(5)*values(5)
                                IF(deriv_idx==1) THEN
                                  WRITE(first_format,"(A,I10,A,I6,A)") "('",nodes_domain%NODES(node_idx)%USER_NUMBER, &
                                    & "',2X,'",deriv_idx,"',5(2X,E12.5))"
                                ELSE
                                  WRITE(first_format,"(A,I6,A)") "(12X,'",deriv_idx,"',5(2X,E12.5))"
                                ENDIF
                                CALL write_string_vector(output_id,1,1,5,5,5,values,first_format,"(20X,5(2X,E12.5))",err,error,*999)
                              ENDDO !deriv_idx
                            ENDDO !node_idx
                            !Output RMS errors                  
                            CALL write_string(output_id,"",err,error,*999)
                            IF(computational_environment%NUMBER_COMPUTATIONAL_NODES>1) THEN
                              IF(any(number>0)) THEN
                                !Local nodes only
                                CALL write_string(output_id,"Local RMS errors:",err,error,*999)
                                local_string= &
                                  & "            Deriv#                                   % error  Absolute err  Relative err"
                                CALL write_string(output_id,local_string,err,error,*999)
                                DO deriv_idx=1,8
                                  IF(number(deriv_idx)>0) THEN
                                    values(1)=sqrt(rms_error_per(deriv_idx)/number(deriv_idx))
                                    values(2)=sqrt(rms_error_abs(deriv_idx)/number(deriv_idx))
                                    values(3)=sqrt(rms_error_rel(deriv_idx)/number(deriv_idx))
                                    WRITE(first_format,"(A,I6,A)") "(12X,'",deriv_idx,"',28X,3(2X,E12.5))"
                                    CALL write_string_vector(output_id,1,1,3,3,3,values,first_format,"(46X,3(2X,E12.5))", &
                                      & err,error,*999)
                                  ENDIF
                                ENDDO !deriv_idx
                                !Local and ghost nodes
                                CALL write_string(output_id,"Local + Ghost RMS errors:",err,error,*999)
                                local_string= &
                                  & "            Deriv#                                   % error  Absolute err  Relative err"
                                CALL write_string(output_id,local_string,err,error,*999)
                                DO deriv_idx=1,8
                                  IF(number(deriv_idx)>0) THEN
                                    values(1)=sqrt((rms_error_per(deriv_idx)+ghost_rms_error_per(deriv_idx))/ &
                                      & (number(deriv_idx)+ghost_number(deriv_idx)))
                                    values(2)=sqrt((rms_error_abs(deriv_idx)+ghost_rms_error_abs(deriv_idx))/ &
                                      & (number(deriv_idx)+ghost_number(deriv_idx)))
                                    values(3)=sqrt((rms_error_rel(deriv_idx)+ghost_rms_error_rel(deriv_idx))/ &
                                      & (number(deriv_idx)+ghost_number(deriv_idx)))
                                    WRITE(first_format,"(A,I6,A)") "(12X,'",deriv_idx,"',28X,3(2X,E12.5))"
                                    CALL write_string_vector(output_id,1,1,3,3,3,values,first_format,"(46X,3(2X,E12.5))", &
                                      & err,error,*999)
                                  ENDIF
                                ENDDO !deriv_idx
                                !Global RMS values
                                !Collect the values across the ranks
                                CALL mpi_allreduce(mpi_in_place,number,8,mpi_integer,mpi_sum,computational_environment%MPI_COMM, &
                                  & mpi_ierror)
                                CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
                                CALL mpi_allreduce(mpi_in_place,rms_error_per,8,mpi_double_precision,mpi_sum, &
                                  & computational_environment%MPI_COMM,mpi_ierror)
                                CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
                                CALL mpi_allreduce(mpi_in_place,rms_error_abs,8,mpi_double_precision,mpi_sum, &
                                  & computational_environment%MPI_COMM,mpi_ierror)
                                CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
                                CALL mpi_allreduce(mpi_in_place,rms_error_rel,8,mpi_double_precision,mpi_sum, &
                                  & computational_environment%MPI_COMM,mpi_ierror)
                                CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
                                CALL write_string(output_id,"Global RMS errors:",err,error,*999)
                                local_string= &
                                  & "            Deriv#                                   % error  Absolute err  Relative err"
                                CALL write_string(output_id,local_string,err,error,*999)
                                DO deriv_idx=1,8
                                  IF(number(deriv_idx)>0) THEN
                                    values(1)=sqrt(rms_error_per(deriv_idx)/number(deriv_idx))
                                    values(2)=sqrt(rms_error_abs(deriv_idx)/number(deriv_idx))
                                    values(3)=sqrt(rms_error_rel(deriv_idx)/number(deriv_idx))
                                    WRITE(first_format,"(A,I6,A)") "(12X,'",deriv_idx,"',28X,3(2X,E12.5))"
                                    CALL write_string_vector(output_id,1,1,3,3,3,values,first_format,"(46X,3(2X,E12.5))", &
                                      & err,error,*999)
                                  ENDIF
                                ENDDO !deriv_idx
                              ENDIF
                            ELSE
                              IF(any(number>0)) THEN
                                CALL write_string(output_id,"RMS errors:",err,error,*999)
                               local_string= &
                                 & "            Deriv#                                   % error  Absolute err  Relative err"
                                CALL write_string(output_id,local_string,err,error,*999)
                                DO deriv_idx=1,8
                                  IF(number(deriv_idx)>0) THEN
                                    values(1)=sqrt(rms_error_per(deriv_idx)/number(deriv_idx))
                                    values(2)=sqrt(rms_error_abs(deriv_idx)/number(deriv_idx))
                                    values(3)=sqrt(rms_error_rel(deriv_idx)/number(deriv_idx))
                                    WRITE(first_format,"(A,I6,A)") "(12X,'",deriv_idx,"',28X,3(2X,E12.5))"
                                    CALL write_string_vector(output_id,1,1,3,3,3,values,first_format,"(46X,3(2X,E12.5))", &
                                      & err,error,*999)
                                  ENDIF
                                ENDDO !deriv_idx
                              ENDIF
                            ENDIF
                          ELSE
                            CALL flagerror("Nodes domain topology is not associated.",err,error,*999)
                          ENDIF                     
                        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
                          local_error="The interpolation type of "// &
                            & trim(number_to_vstring(field%VARIABLES(var_idx)%COMPONENTS(component_idx)% &
                            & interpolation_type,"*",err,error))//" for component number "// &
                            & trim(number_to_vstring(component_idx,"*",err,error))//" of variable type "// &
                            & trim(number_to_vstring(variable_type,"*",err,error))//" of field number "// &
                            & trim(number_to_vstring(field%USER_NUMBER,"*",err,error))//" is invalid."
                          CALL flagerror(local_error,err,error,*999)
                        END SELECT
                      ELSE
                        CALL flagerror("Domain topology is not associated.",err,error,*999)
                      ENDIF
                    ELSE
                      CALL flagerror("Domain is not associated.",err,error,*999)
                    ENDIF
                    CALL write_string(output_id,"",err,error,*999)
                  ENDDO !component_idx
                  !Restore the dependent and analytic parameter sets
                  CALL field_parameter_set_data_restore(field,variable_type,field_values_set_type,numerical_values, &
                    & err,error,*999)
                  CALL field_parameter_set_data_restore(field,variable_type,field_analytic_values_set_type,analytic_values, &
                    & err,error,*999)
                  !Allocated the integral errors
                  ALLOCATE(integral_errors(6,field_variable%NUMBER_OF_COMPONENTS),stat=err)
                  IF(err/=0) CALL flagerror("Could not allocate integral errors.",err,error,*999)
                  ALLOCATE(ghost_integral_errors(6,field_variable%NUMBER_OF_COMPONENTS),stat=err)
                  IF(err/=0) CALL flagerror("Could not allocate ghost integral errors.",err,error,*999)
                  CALL analytic_analysis_integral_errors(field_variable,integral_errors,ghost_integral_errors,err,error,*999)
                  IF(computational_environment%NUMBER_COMPUTATIONAL_NODES>1) THEN
                    CALL write_string(output_id,"Local Integral errors:",err,error,*999)
                    local_string="Component#             Numerical      Analytic       % error  Absolute err  Relative err"
                    CALL write_string(output_id,local_string,err,error,*999)
                    DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
                      values(1)=integral_errors(1,component_idx)
                      values(2)=integral_errors(3,component_idx)
                      values(3)=analytic_analysis_percentage_error(values(1),values(2))
                      values(4)=analytic_analysis_absolute_error(values(1),values(2))
                      values(5)=analytic_analysis_relative_error(values(1),values(2))
                      WRITE(first_format,"(A,I10,A)") "('",component_idx,"',2X,'Intg  ',5(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,5,5,5,values,first_format,"(20X,5(2X,E12.5))",err,error,*999)
                      values(1)=integral_errors(2,component_idx)
                      values(2)=integral_errors(4,component_idx)
                      values(3)=analytic_analysis_percentage_error(values(1),values(2))
                      values(4)=analytic_analysis_absolute_error(values(1),values(2))
                      values(5)=analytic_analysis_relative_error(values(1),values(2))
                      WRITE(first_format,"(A)") "(12X,'Intg^2',5(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,5,5,5,values,first_format,"(20X,5(2X,E12.5))",err,error,*999)
                    ENDDO !component_idx
                    local_string="Component#             Numerical      Analytic           NID        NID(%)"
                    CALL write_string(output_id,local_string,err,error,*999)
                    DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
                      values(1)=integral_errors(5,component_idx)
                      values(2)=integral_errors(3,component_idx)
                      values(3)=analytic_analysis_nid_error(values(1),values(2))
                      values(4)=values(3)*100.0_dp
                     WRITE(first_format,"(A,I10,A)") "('",component_idx,"',2X,'Diff  ',4(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,4,4,4,values,first_format,"(20X,4(2X,E12.5))",err,error,*999)
                      values(1)=integral_errors(6,component_idx)
                      values(2)=integral_errors(4,component_idx)
                      values(3)=analytic_analysis_nid_error(values(1),values(2))
                      values(4)=values(3)*100.0_dp
                      WRITE(first_format,"(A)") "(12X,'Diff^2',4(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,4,4,4,values,first_format,"(20X,4(2X,E12.5))",err,error,*999)
                    ENDDO !component_idx
                    CALL write_string(output_id,"Local + Ghost Integral errors:",err,error,*999)
                    local_string="Component#             Numerical      Analytic       % error  Absolute err  Relative err"
                    CALL write_string(output_id,local_string,err,error,*999)
                    DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
                      values(1)=integral_errors(1,component_idx)+ghost_integral_errors(1,component_idx)
                      values(2)=integral_errors(3,component_idx)+ghost_integral_errors(3,component_idx)
                      values(3)=analytic_analysis_percentage_error(values(1),values(2))
                      values(4)=analytic_analysis_absolute_error(values(1),values(2))
                      values(5)=analytic_analysis_relative_error(values(1),values(2))
                      WRITE(first_format,"(A,I10,A)") "('",component_idx,"',2X,'Intg  ',5(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,5,5,5,values,first_format,"(20X,5(2X,E12.5))",err,error,*999)
                      values(1)=integral_errors(2,component_idx)
                      values(2)=integral_errors(4,component_idx)
                      values(3)=analytic_analysis_percentage_error(values(1),values(2))
                      values(4)=analytic_analysis_absolute_error(values(1),values(2))
                      values(5)=analytic_analysis_relative_error(values(1),values(2))
                      WRITE(first_format,"(A)") "(12X,'Intg^2',5(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,5,5,5,values,first_format,"(20X,5(2X,E12.5))",err,error,*999)
                    ENDDO !component_idx
                    local_string="Component#             Numerical      Analytic           NID        NID(%)"
                    CALL write_string(output_id,local_string,err,error,*999)
                    DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
                      values(1)=integral_errors(5,component_idx)+ghost_integral_errors(5,component_idx)
                      values(2)=integral_errors(3,component_idx)+ghost_integral_errors(3,component_idx)
                      values(3)=analytic_analysis_nid_error(values(1),values(2))
                      values(4)=values(3)*100.0_dp
                      WRITE(first_format,"(A,I10,A)") "('",component_idx,"',2X,'Diff  ',4(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,4,4,4,values,first_format,"(20X,4(2X,E12.5))",err,error,*999)
                      values(1)=integral_errors(6,component_idx)+ghost_integral_errors(6,component_idx)
                      values(2)=integral_errors(4,component_idx)+ghost_integral_errors(4,component_idx)
                      values(3)=analytic_analysis_nid_error(values(1),values(2))
                      values(4)=values(3)*100.0_dp
                      WRITE(first_format,"(A)") "(12X,'Diff^2',4(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,4,4,4,values,first_format,"(20X,4(2X,E12.5))",err,error,*999)
                    ENDDO !component_idx
                    !Collect the values across the ranks
                    CALL mpi_allreduce(mpi_in_place,integral_errors,6*field_variable%NUMBER_OF_COMPONENTS,mpi_double_precision, &
                      & mpi_sum,computational_environment%MPI_COMM,mpi_ierror)
                    CALL write_string(output_id,"Global Integral errors:",err,error,*999)
                    local_string="Component#             Numerical      Analytic       % error  Absolute err  Relative err"
                    CALL write_string(output_id,local_string,err,error,*999)
                    DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
                      CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
                      values(1)=integral_errors(1,component_idx)
                      values(2)=integral_errors(3,component_idx)
                      values(4)=analytic_analysis_absolute_error(values(1),values(2))
                      values(5)=analytic_analysis_relative_error(values(1),values(2))
                      WRITE(first_format,"(A,I10,A)") "('",component_idx,"',2X,'Intg  ',5(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,5,5,5,values,first_format,"(20X,5(2X,E12.5))",err,error,*999)
                      values(1)=integral_errors(2,component_idx)
                      values(2)=integral_errors(4,component_idx)
                      values(3)=analytic_analysis_percentage_error(values(1),values(2))
                      values(4)=analytic_analysis_absolute_error(values(1),values(2))
                      values(5)=analytic_analysis_relative_error(values(1),values(2))
                      WRITE(first_format,"(A)") "(12X,'Intg^2',5(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,5,5,5,values,first_format,"(20X,5(2X,E12.5))",err,error,*999) 
                    ENDDO !component_idx
                    local_string="Component#             Numerical      Analytic           NID        NID(%)"
                    CALL write_string(output_id,local_string,err,error,*999)
                    DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
                      values(1)=integral_errors(5,component_idx)
                      values(2)=integral_errors(3,component_idx)
                      values(3)=analytic_analysis_nid_error(values(1),values(2))
                      values(4)=values(3)*100.0_dp
                     WRITE(first_format,"(A,I10,A)") "('",component_idx,"',2X,'Diff  ',4(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,4,4,4,values,first_format,"(20X,4(2X,E12.5))",err,error,*999)
                      values(1)=integral_errors(6,component_idx)
                      values(2)=integral_errors(4,component_idx)
                      values(3)=analytic_analysis_nid_error(values(1),values(2))
                      values(4)=values(3)*100.0_dp
                      WRITE(first_format,"(A)") "(12X,'Diff^2',4(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,4,4,4,values,first_format,"(20X,4(2X,E12.5))",err,error,*999)
                    ENDDO !component_idx
                  ELSE
                    CALL write_string(output_id,"Integral errors:",err,error,*999)
                    local_string="Component#             Numerical      Analytic       % error  Absolute err  Relative err"
                    CALL write_string(output_id,local_string,err,error,*999)
                    DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
                      values(1)=integral_errors(1,component_idx)
                      values(2)=integral_errors(3,component_idx)
                      values(3)=analytic_analysis_percentage_error(values(1),values(2))
                      values(4)=analytic_analysis_absolute_error(values(1),values(2))
                      values(5)=analytic_analysis_relative_error(values(1),values(2))
                      WRITE(first_format,"(A,I10,A)") "('",component_idx,"',2X,'Intg  ',5(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,5,5,5,values,first_format,"(20X,5(2X,E12.5))",err,error,*999)
                      values(1)=integral_errors(2,component_idx)
                      values(2)=integral_errors(4,component_idx)
                      values(3)=analytic_analysis_percentage_error(values(1),values(2))
                      values(4)=analytic_analysis_absolute_error(values(1),values(2))
                      values(5)=analytic_analysis_relative_error(values(1),values(2))
                      WRITE(first_format,"(A)") "(12X,'Intg^2',5(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,5,5,5,values,first_format,"(20X,5(2X,E12.5))",err,error,*999)
                    ENDDO !component_idx
                    local_string="Component#             Numerical      Analytic           NID        NID(%)"
                    CALL write_string(output_id,local_string,err,error,*999)
                    DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
                      values(1)=integral_errors(5,component_idx)
                      values(2)=integral_errors(3,component_idx)
                      values(3)=analytic_analysis_nid_error(values(1),values(2))
                      values(4)=values(3)*100.0_dp
                     WRITE(first_format,"(A,I10,A)") "('",component_idx,"',2X,'Diff  ',4(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,4,4,4,values,first_format,"(20X,4(2X,E12.5))",err,error,*999)
                      values(1)=integral_errors(6,component_idx)
                      values(2)=integral_errors(4,component_idx)
                      values(3)=analytic_analysis_nid_error(values(1),values(2))
                      values(4)=values(3)*100.0_dp
                      WRITE(first_format,"(A)") "(12X,'Diff^2',4(2X,E12.5))"
                      CALL write_string_vector(output_id,1,1,4,4,4,values,first_format,"(20X,4(2X,E12.5))",err,error,*999)
                    ENDDO !component_idx
                    CALL write_string(output_id,"",err,error,*999)
                  ENDIF
                  IF(ALLOCATED(integral_errors)) DEALLOCATE(integral_errors)
                  IF(ALLOCATED(ghost_integral_errors)) DEALLOCATE(ghost_integral_errors)
                ELSE
                  CALL flagerror("Field variable is not associated.",err,error,*999)
                ENDIF
              ENDDO !var_idx
            ELSE
              CALL flagerror("Decomposition topology is not associated.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Field decomposition is not associated.",err,error,*999)
          ENDIF
          IF(len_trim(filename)>=1) THEN
            CLOSE(unit=output_id,iostat=err)
            IF(err/=0) CALL flagerror("Error closing analysis output file.",err,error,*999)
          ENDIF
        ELSE
          local_error="Field number "//trim(number_to_vstring(field%USER_NUMBER,"*",err,error))// &
            & " is not a dependent field."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        local_error="Field number "//trim(number_to_vstring(field%USER_NUMBER,"*",err,error))// &
          & " has not been finished."
        CALL flagerror(local_error,err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Field is not associated.",err,error,*999)
    ENDIF
          
    exits("AnalyticAnalysis_Output")
    RETURN
999 IF(ALLOCATED(integral_errors)) DEALLOCATE(integral_errors)
    IF(ALLOCATED(ghost_integral_errors)) DEALLOCATE(ghost_integral_errors)
    errorsexits("AnalyticAnalysis_Output",err,error)
    RETURN 1
  END SUBROUTINE analyticanalysis_output

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

  PURE FUNCTION analytic_analysis_absolute_error(NUMERIC_VALUE,ANALYTIC_VALUE)

    !Argument variables
    REAL(DP), INTENT(IN) :: NUMERIC_VALUE
    REAL(DP), INTENT(IN) :: ANALYTIC_VALUE
    !Function result
    REAL(DP) :: ANALYTIC_ANALYSIS_ABSOLUTE_ERROR

    analytic_analysis_absolute_error=abs(analytic_value-numeric_value)

  END FUNCTION analytic_analysis_absolute_error
  
  !
  !================================================================================================================================
  !  

  PURE FUNCTION analytic_analysis_percentage_error(NUMERIC_VALUE,ANALYTIC_VALUE)

    !Argument variables
    REAL(DP), INTENT(IN) :: NUMERIC_VALUE
    REAL(DP), INTENT(IN) :: ANALYTIC_VALUE
    !Function result
    REAL(DP) :: ANALYTIC_ANALYSIS_PERCENTAGE_ERROR

    IF(abs(analytic_value)>zero_tolerance) THEN
      analytic_analysis_percentage_error=(analytic_value-numeric_value)/analytic_value*100.0_dp
    ELSE
      analytic_analysis_percentage_error=0.0_dp
    ENDIF

  END FUNCTION analytic_analysis_percentage_error
  
  !
  !================================================================================================================================
  !  

  PURE FUNCTION analytic_analysis_relative_error(NUMERIC_VALUE,ANALYTIC_VALUE)

    !Argument variables
    REAL(DP), INTENT(IN) :: NUMERIC_VALUE
    REAL(DP), INTENT(IN) :: ANALYTIC_VALUE
    !Function result
    REAL(DP) :: ANALYTIC_ANALYSIS_RELATIVE_ERROR

    IF(abs(1.0_dp+analytic_value)>zero_tolerance) THEN
      analytic_analysis_relative_error=(analytic_value-numeric_value)/(1.0_dp+analytic_value)
    ELSE
      analytic_analysis_relative_error=0.0_dp
    ENDIF

  END FUNCTION analytic_analysis_relative_error
  
  !
  !================================================================================================================================
  !  

  PURE FUNCTION analytic_analysis_nid_error(NUMERIC_VALUE,ANALYTIC_VALUE)

    !Argument variables
    REAL(DP), INTENT(IN) :: NUMERIC_VALUE
    REAL(DP), INTENT(IN) :: ANALYTIC_VALUE
    !Function result
    REAL(DP) :: ANALYTIC_ANALYSIS_NID_ERROR

    IF(abs(analytic_value)>zero_tolerance) THEN
      analytic_analysis_nid_error=(analytic_value-numeric_value)/analytic_value
    ELSE
      analytic_analysis_nid_error=(analytic_value-numeric_value)/(1.0_dp+analytic_value)
    ENDIF

  END FUNCTION analytic_analysis_nid_error
  
  !
  !================================================================================================================================
  !  

  SUBROUTINE analytic_analysis_integral_errors(FIELD_VARIABLE,INTEGRAL_ERRORS,GHOST_INTEGRAL_ERRORS,ERR,ERROR,*)

    !Argument variables
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
    REAL(DP), INTENT(OUT) :: INTEGRAL_ERRORS(:,:)
    REAL(DP), INTENT(OUT) :: GHOST_INTEGRAL_ERRORS(:,:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: element_idx,component_idx,gauss_idx,parameter_idx,variable_type
    REAL(DP) :: ANALYTIC_INT,NUMERICAL_INT,RWG
    TYPE(basis_type), POINTER :: BASIS,DEPENDENT_BASIS,GEOMETRIC_BASIS
    TYPE(decomposition_type), POINTER :: DECOMPOSITION
    TYPE(domain_elements_type), POINTER :: DOMAIN_ELEMENTS1,DOMAIN_ELEMENTS2,DOMAIN_ELEMENTS3
    TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD
    TYPE(field_interpolated_point_ptr_type), POINTER :: GEOMETRIC_INTERP_POINT(:)
    TYPE(field_interpolated_point_metrics_ptr_type), POINTER :: GEOMETRIC_INTERP_POINT_METRICS(:)
    TYPE(field_interpolation_parameters_ptr_type), POINTER :: ANALYTIC_INTERP_PARAMETERS(:),GEOMETRIC_INTERP_PARAMETERS(:), &
      & NUMERICAL_INTERP_PARAMETERS(:)
    TYPE(field_variable_type), POINTER :: GEOMETRIC_VARIABLE
    TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME
    TYPE(varying_string) :: LOCAL_ERROR

    NULLIFY(geometric_interp_point)
    NULLIFY(geometric_interp_point_metrics)
    NULLIFY(analytic_interp_parameters)
    NULLIFY(geometric_interp_parameters)
    NULLIFY(numerical_interp_parameters)

    enters("ANALYTIC_ANALYSIS_INTEGRAL_ERRORS",err,error,*999)
    
    integral_errors=0.0_dp
    ghost_integral_errors=0.0_dp
    IF(ASSOCIATED(field_variable)) THEN
      IF(SIZE(integral_errors,1)>=6.AND.SIZE(integral_errors,2)>=field_variable%NUMBER_OF_COMPONENTS) THEN
        IF(SIZE(ghost_integral_errors,1)>=6.AND.SIZE(ghost_integral_errors,2)>=field_variable%NUMBER_OF_COMPONENTS) THEN
          variable_type=field_variable%VARIABLE_TYPE
          dependent_field=>field_variable%FIELD
          IF(ASSOCIATED(dependent_field)) THEN
            decomposition=>dependent_field%DECOMPOSITION
            IF(ASSOCIATED(decomposition)) THEN
              geometric_field=>dependent_field%GEOMETRIC_FIELD
              IF(ASSOCIATED(geometric_field)) THEN
                geometric_variable=>geometric_field%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR
                IF(ASSOCIATED(geometric_variable)) THEN
                  CALL field_interpolation_parameters_initialise(geometric_field,geometric_interp_parameters,err,error,*999)
                  CALL field_interpolation_parameters_initialise(dependent_field,numerical_interp_parameters,err,error,*999)
                  CALL field_interpolation_parameters_initialise(dependent_field,analytic_interp_parameters,err,error,*999)
                  CALL field_interpolated_points_initialise(geometric_interp_parameters,geometric_interp_point,err,error,*999)
                  CALL field_interpolatedpointsmetricsinitialise(geometric_interp_point,geometric_interp_point_metrics, &
                    & err,error,*999)
                  domain_elements1=>field_variable%COMPONENTS(decomposition%MESH_COMPONENT_NUMBER)%DOMAIN%TOPOLOGY%ELEMENTS
                  domain_elements2=>geometric_variable%COMPONENTS(decomposition%MESH_COMPONENT_NUMBER)%DOMAIN%TOPOLOGY%ELEMENTS
                  DO element_idx=1,domain_elements1%NUMBER_OF_ELEMENTS
                    dependent_basis=>domain_elements1%ELEMENTS(element_idx)%BASIS
                    geometric_basis=>domain_elements2%ELEMENTS(element_idx)%BASIS
                    quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
                    CALL field_interpolation_parameters_element_get(field_values_set_type,element_idx, &
                      & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
                    CALL field_interpolation_parameters_element_get(field_values_set_type,element_idx, &
                      & numerical_interp_parameters(variable_type)%PTR,err,error,*999)
                    CALL field_interpolation_parameters_element_get(field_analytic_values_set_type,element_idx, &
                      & analytic_interp_parameters(variable_type)%PTR,err,error,*999)
                    DO gauss_idx=1,quadrature_scheme%NUMBER_OF_GAUSS
                      CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,gauss_idx, &
                        & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
                      CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI, &
                        & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
                      rwg=geometric_interp_point_metrics(field_u_variable_type)%PTR%JACOBIAN* &
                        & quadrature_scheme%GAUSS_WEIGHTS(gauss_idx)
                      DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
                        domain_elements3=>field_variable%COMPONENTS(component_idx)%DOMAIN%TOPOLOGY%ELEMENTS
                        basis=>domain_elements3%ELEMENTS(element_idx)%BASIS
                        numerical_int=0.0_dp
                        analytic_int=0.0_dp
                        SELECT CASE(dependent_field%SCALINGS%SCALING_TYPE)
                        CASE(field_no_scaling)
                          DO parameter_idx=1,basis%NUMBER_OF_ELEMENT_PARAMETERS
                            numerical_int=numerical_int+quadrature_scheme%GAUSS_BASIS_FNS(parameter_idx,no_part_deriv,gauss_idx)* &
                              & numerical_interp_parameters(variable_type)%PTR%PARAMETERS(parameter_idx,component_idx)
                            analytic_int=analytic_int+quadrature_scheme%GAUSS_BASIS_FNS(parameter_idx,no_part_deriv,gauss_idx)* &
                              & analytic_interp_parameters(variable_type)%PTR%PARAMETERS(parameter_idx,component_idx)
                          ENDDO !parameter_idx
                        CASE(field_unit_scaling,field_arc_length_scaling,field_arithmetic_mean_scaling,field_harmonic_mean_scaling)
                          DO parameter_idx=1,basis%NUMBER_OF_ELEMENT_PARAMETERS
                            numerical_int=numerical_int+quadrature_scheme%GAUSS_BASIS_FNS(parameter_idx,no_part_deriv,gauss_idx)* &
                              & numerical_interp_parameters(variable_type)%PTR%PARAMETERS(parameter_idx,component_idx)* &
                              & numerical_interp_parameters(variable_type)%PTR%SCALE_FACTORS(parameter_idx,component_idx)
                            analytic_int=analytic_int+quadrature_scheme%GAUSS_BASIS_FNS(parameter_idx,no_part_deriv,gauss_idx)* &
                              & analytic_interp_parameters(variable_type)%PTR%PARAMETERS(parameter_idx,component_idx)* &
                              & analytic_interp_parameters(variable_type)%PTR%SCALE_FACTORS(parameter_idx,component_idx)
                          ENDDO !parameter_idx
                        CASE DEFAULT
                          local_error="The dependent field scaling type of "// &
                            & trim(number_to_vstring(dependent_field%SCALINGS%SCALING_TYPE,"*",err,error))//" is invalid."
                          CALL flagerror(local_error,err,error,*999)
                        END SELECT
                        integral_errors(1,component_idx)=integral_errors(1,component_idx)+numerical_int*rwg
                        integral_errors(2,component_idx)=integral_errors(2,component_idx)+numerical_int**2*rwg
                        integral_errors(3,component_idx)=integral_errors(3,component_idx)+analytic_int*rwg
                        integral_errors(4,component_idx)=integral_errors(4,component_idx)+analytic_int**2*rwg
                        integral_errors(5,component_idx)=integral_errors(5,component_idx)+(analytic_int-numerical_int)*rwg
                        integral_errors(6,component_idx)=integral_errors(6,component_idx)+(analytic_int-numerical_int)**2*rwg
                      ENDDO !component_idx
                    ENDDO !gauss_idx
                  ENDDO !element_idx
                  DO element_idx=domain_elements1%NUMBER_OF_ELEMENTS+1,domain_elements1%TOTAL_NUMBER_OF_ELEMENTS
                    dependent_basis=>domain_elements1%ELEMENTS(element_idx)%BASIS
                    geometric_basis=>domain_elements2%ELEMENTS(element_idx)%BASIS
                    quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
                    CALL field_interpolation_parameters_element_get(field_values_set_type,element_idx, &
                      & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
                    CALL field_interpolation_parameters_element_get(field_values_set_type,element_idx, &
                      & numerical_interp_parameters(variable_type)%PTR,err,error,*999)
                    CALL field_interpolation_parameters_element_get(field_analytic_values_set_type,element_idx, &
                      & analytic_interp_parameters(variable_type)%PTR,err,error,*999)
                    DO gauss_idx=1,quadrature_scheme%NUMBER_OF_GAUSS
                      CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,gauss_idx, &
                        & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
                      CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI, &
                        & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
                      rwg=geometric_interp_point_metrics(field_u_variable_type)%PTR%JACOBIAN* &
                        & quadrature_scheme%GAUSS_WEIGHTS(gauss_idx)
                      DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
                        domain_elements3=>field_variable%COMPONENTS(component_idx)%DOMAIN%TOPOLOGY%ELEMENTS
                        basis=>domain_elements3%ELEMENTS(element_idx)%BASIS
                        numerical_int=0.0_dp
                        analytic_int=0.0_dp
                        SELECT CASE(dependent_field%SCALINGS%SCALING_TYPE)
                        CASE(field_no_scaling)
                          DO parameter_idx=1,basis%NUMBER_OF_ELEMENT_PARAMETERS
                            numerical_int=numerical_int+quadrature_scheme%GAUSS_BASIS_FNS(parameter_idx,no_part_deriv,gauss_idx)* &
                              & numerical_interp_parameters(variable_type)%PTR%PARAMETERS(parameter_idx,component_idx)
                            analytic_int=analytic_int+quadrature_scheme%GAUSS_BASIS_FNS(parameter_idx,no_part_deriv,gauss_idx)* &
                              & analytic_interp_parameters(variable_type)%PTR%PARAMETERS(parameter_idx,component_idx)
                          ENDDO !parameter_idx
                        CASE(field_unit_scaling,field_arc_length_scaling,field_arithmetic_mean_scaling,field_harmonic_mean_scaling)
                          DO parameter_idx=1,basis%NUMBER_OF_ELEMENT_PARAMETERS
                            numerical_int=numerical_int+quadrature_scheme%GAUSS_BASIS_FNS(parameter_idx,no_part_deriv,gauss_idx)* &
                              & numerical_interp_parameters(variable_type)%PTR%PARAMETERS(parameter_idx,component_idx)* &
                              & numerical_interp_parameters(variable_type)%PTR%SCALE_FACTORS(parameter_idx,component_idx)
                            analytic_int=analytic_int+quadrature_scheme%GAUSS_BASIS_FNS(parameter_idx,no_part_deriv,gauss_idx)* &
                              & analytic_interp_parameters(variable_type)%PTR%PARAMETERS(parameter_idx,component_idx)* &
                              & analytic_interp_parameters(variable_type)%PTR%SCALE_FACTORS(parameter_idx,component_idx)
                          ENDDO !parameter_idx
                        CASE DEFAULT
                          local_error="The dependent field scaling type of "// &
                            & trim(number_to_vstring(dependent_field%SCALINGS%SCALING_TYPE,"*",err,error))//" is invalid."
                          CALL flagerror(local_error,err,error,*999)
                        END SELECT
                        ghost_integral_errors(1,component_idx)=ghost_integral_errors(1,component_idx)+numerical_int*rwg
                        ghost_integral_errors(2,component_idx)=ghost_integral_errors(2,component_idx)+numerical_int**2*rwg
                        ghost_integral_errors(3,component_idx)=ghost_integral_errors(3,component_idx)+analytic_int*rwg
                        ghost_integral_errors(4,component_idx)=ghost_integral_errors(4,component_idx)+analytic_int**2*rwg
                        ghost_integral_errors(5,component_idx)=ghost_integral_errors(5,component_idx)+ &
                          & (analytic_int-numerical_int)*rwg
                        ghost_integral_errors(6,component_idx)=ghost_integral_errors(6,component_idx)+ &
                          & (analytic_int-numerical_int)**2*rwg
                      ENDDO !component_idx
                    ENDDO !gauss_idx
                  ENDDO !element_idx
                  CALL field_interpolatedpointsmetricsfinalise(geometric_interp_point_metrics,err,error,*999)
                  CALL field_interpolated_points_finalise(geometric_interp_point,err,error,*999)
                  CALL field_interpolation_parameters_finalise(analytic_interp_parameters,err,error,*999)
                  CALL field_interpolation_parameters_finalise(numerical_interp_parameters,err,error,*999)
                  CALL field_interpolation_parameters_finalise(geometric_interp_parameters,err,error,*999)
                ELSE
                  CALL flagerror("Geometric field variable is not associated.",err,error,*999)
                ENDIF
              ELSE
                CALL flagerror("Field geometric field is not associated.",err,error,*999)
              ENDIF
            ELSE
              CALL flagerror("Field decomposition is not associated.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Field variable field is not associated.",err,error,*999)
          ENDIF
        ELSE
          local_error="Invalid size for GHOST_INTEGRAL_ERRORS. The size is ("// &
            & trim(number_to_vstring(SIZE(ghost_integral_errors,1),"*",err,error))//","// &
            & trim(number_to_vstring(SIZE(ghost_integral_errors,2),"*",err,error))//") and it needs to be at least (6,"// &
            & trim(number_to_vstring(field_variable%NUMBER_OF_COMPONENTS,"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        local_error="Invalid size for INTEGRAL_ERRORS. The size is ("// &
          & trim(number_to_vstring(SIZE(integral_errors,1),"*",err,error))//","// &
          & trim(number_to_vstring(SIZE(integral_errors,2),"*",err,error))//") and it needs to be at least (6,"// &
          & trim(number_to_vstring(field_variable%NUMBER_OF_COMPONENTS,"*",err,error))//")."
        CALL flagerror(local_error,err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Field variable is not associated.",err,error,*999)
    ENDIF

    exits("ANALYTIC_ANALYSIS_INTEGRAL_ERRORS")
    RETURN
999 errorsexits("ANALYTIC_ANALYSIS_INTEGRAL_ERRORS",err,error)
    RETURN 1
  END SUBROUTINE analytic_analysis_integral_errors

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

  SUBROUTINE analyticanalysis_integralabsoluteerrorget(FIELD,VARIABLE_TYPE,COMPONENT_NUMBER,INTEGRAL_ERROR, &
    & ghost_integral_error,err,error,*)

    !Argument variables   
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    REAL(DP), INTENT(OUT) :: INTEGRAL_ERROR(2)
    REAL(DP), INTENT(OUT) :: GHOST_INTEGRAL_ERROR(2)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
    REAL(DP), ALLOCATABLE :: INTEGRAL_ERRORS(:,:)
    REAL(DP), ALLOCATABLE :: GHOST_INTEGRAL_ERRORS(:,:)

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

    IF(ASSOCIATED(field)) THEN
      field_variable=>field%VARIABLE_TYPE_MAP(variable_type)%PTR
      CALL analytic_analysis_integral_errors(field_variable,integral_errors,ghost_integral_errors,err,error,*999)
      integral_error(1)=analytic_analysis_absolute_error(integral_errors(1,component_number),integral_errors(3,component_number))
      integral_error(2)=analytic_analysis_absolute_error(integral_errors(2,component_number),integral_errors(4,component_number))
      ghost_integral_error(1)=analytic_analysis_absolute_error(ghost_integral_errors(1,component_number), &
        & ghost_integral_errors(3,component_number))
      ghost_integral_error(2)=analytic_analysis_absolute_error(ghost_integral_errors(2,component_number), &
        & ghost_integral_errors(4,component_number))
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF 
    
    exits("AnalyticAnalysis_IntegralAbsoluteErrorGet")
    RETURN
999 errorsexits("AnalyticAnalysis_IntegralAbsoluteErrorGet",err,error)
    RETURN 1
    
  END SUBROUTINE analyticanalysis_integralabsoluteerrorget
  
   !
  !================================================================================================================================
  !

  SUBROUTINE analyticanalysis_integralanalyticvalueget(FIELD,VARIABLE_TYPE,COMPONENT_NUMBER,INTEGRAL_ERROR, &
    & ghost_integral_error,err,error,*)
  
    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    REAL(DP), INTENT(OUT) :: INTEGRAL_ERROR(2)
    REAL(DP), INTENT(OUT) :: GHOST_INTEGRAL_ERROR(2)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
    REAL(DP), ALLOCATABLE :: INTEGRAL_ERRORS(:,:)
    REAL(DP), ALLOCATABLE :: GHOST_INTEGRAL_ERRORS(:,:)

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

    IF(ASSOCIATED(field)) THEN
      field_variable=>field%VARIABLE_TYPE_MAP(variable_type)%PTR
      CALL analytic_analysis_integral_errors(field_variable,integral_errors,ghost_integral_errors,err,error,*999)
      integral_error(1)=integral_errors(3,component_number)
      integral_error(2)=integral_errors(4,component_number)
      ghost_integral_error(1)=ghost_integral_errors(3,component_number)
      ghost_integral_error(2)=ghost_integral_errors(4,component_number)
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF 
    
    exits("AnalyticAnalysis_IntegralAnalyticValueGet")
    RETURN
999 errorsexits("AnalyticAnalysis_IntegralAnalyticValueGet",err,error)
    RETURN 1
    
  END SUBROUTINE analyticanalysis_integralanalyticvalueget
  
  !
  !================================================================================================================================
  !

  SUBROUTINE analyticanalysis_integralnumericalvalueget(FIELD,VARIABLE_TYPE,COMPONENT_NUMBER,INTEGRAL_ERROR, &
    & ghost_integral_error,err,error,*)
  
    !Argument variables   
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    REAL(DP), INTENT(OUT) :: INTEGRAL_ERROR(2)
    REAL(DP), INTENT(OUT) :: GHOST_INTEGRAL_ERROR(2)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
    REAL(DP), ALLOCATABLE :: INTEGRAL_ERRORS(:,:)
    REAL(DP), ALLOCATABLE :: GHOST_INTEGRAL_ERRORS(:,:)

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

    IF(ASSOCIATED(field)) THEN
      field_variable=>field%VARIABLE_TYPE_MAP(variable_type)%PTR
      CALL analytic_analysis_integral_errors(field_variable,integral_errors,ghost_integral_errors,err,error,*999)
      integral_error(1)=integral_errors(1,component_number)
      integral_error(2)=integral_errors(2,component_number)
      ghost_integral_error(1)=ghost_integral_errors(1,component_number)
      ghost_integral_error(2)=ghost_integral_errors(2,component_number)
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF 
    
    exits("AnalyticAnalysis_IntegralNumericalValueGet")
    RETURN
999 errorsexits("AnalyticAnalysis_IntegralNumericalValueGet",err,error)
    RETURN 1
    
  END SUBROUTINE analyticanalysis_integralnumericalvalueget
  
  !
  !================================================================================================================================
  !

  SUBROUTINE analyticanalysis_integralnidnumericalvalueget(FIELD,VARIABLE_TYPE,COMPONENT_NUMBER,INTEGRAL_ERROR, &
    & ghost_integral_error,err,error,*)

    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    REAL(DP), INTENT(OUT) :: INTEGRAL_ERROR(2)
    REAL(DP), INTENT(OUT) :: GHOST_INTEGRAL_ERROR(2)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
    REAL(DP), ALLOCATABLE :: INTEGRAL_ERRORS(:,:)
    REAL(DP), ALLOCATABLE :: GHOST_INTEGRAL_ERRORS(:,:)

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

    IF(ASSOCIATED(field)) THEN
      field_variable=>field%VARIABLE_TYPE_MAP(variable_type)%PTR
      CALL analytic_analysis_integral_errors(field_variable,integral_errors,ghost_integral_errors,err,error,*999)
      integral_error(1)=integral_errors(5,component_number)
      integral_error(2)=integral_errors(6,component_number)
      ghost_integral_error(1)=ghost_integral_errors(5,component_number)
      ghost_integral_error(2)=ghost_integral_errors(6,component_number)
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF

    exits("AnalyticAnalysis_IntegralNIDNumericalValueGet")
    RETURN
999 errors("AnalyticAnalysis_IntegralNIDNumericalValueGet",err,error)
    exits("AnalyticAnalysis_IntegralNIDNumericalValueGet")
    RETURN 1
    
  END SUBROUTINE analyticanalysis_integralnidnumericalvalueget

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

  SUBROUTINE analyticanalysis_integralniderrorget(FIELD,VARIABLE_TYPE,COMPONENT_NUMBER,INTEGRAL_ERROR, &
    & ghost_integral_error,err,error,*)

    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    REAL(DP), INTENT(OUT) :: INTEGRAL_ERROR(2)
    REAL(DP), INTENT(OUT) :: GHOST_INTEGRAL_ERROR(2)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
    REAL(DP), ALLOCATABLE :: INTEGRAL_ERRORS(:,:)
    REAL(DP), ALLOCATABLE :: GHOST_INTEGRAL_ERRORS(:,:)

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

    IF(ASSOCIATED(field)) THEN
      field_variable=>field%VARIABLE_TYPE_MAP(variable_type)%PTR
      CALL analytic_analysis_integral_errors(field_variable,integral_errors,ghost_integral_errors,err,error,*999)
      integral_error(1)=analytic_analysis_nid_error(integral_errors(5,component_number),integral_errors(3,component_number))
      integral_error(2)=analytic_analysis_nid_error(integral_errors(6,component_number),integral_errors(4,component_number))
      ghost_integral_error(1)=analytic_analysis_nid_error(ghost_integral_errors(5,component_number), &
        & ghost_integral_errors(3,component_number))
      ghost_integral_error(2)=analytic_analysis_nid_error(ghost_integral_errors(6,component_number), &
        & ghost_integral_errors(4,component_number))
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF

    exits("AnalyticAnalysis_IntegralNIDErrorGet")
    RETURN
999 errorsexits("AnalyticAnalysis_IntegralNIDErrorGet",err,error)
    RETURN 1
  END SUBROUTINE analyticanalysis_integralniderrorget

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

  SUBROUTINE analyticanalysis_integralpercentageerrorget(FIELD,VARIABLE_TYPE,COMPONENT_NUMBER,INTEGRAL_ERROR, &
    & ghost_integral_error,err,error,*)

    !Argument variables   
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    REAL(DP), INTENT(OUT) :: INTEGRAL_ERROR(2)
    REAL(DP), INTENT(OUT) :: GHOST_INTEGRAL_ERROR(2)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
    REAL(DP), ALLOCATABLE :: INTEGRAL_ERRORS(:,:)
    REAL(DP), ALLOCATABLE :: GHOST_INTEGRAL_ERRORS(:,:)

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

    IF(ASSOCIATED(field)) THEN
      field_variable=>field%VARIABLE_TYPE_MAP(variable_type)%PTR
      CALL analytic_analysis_integral_errors(field_variable,integral_errors,ghost_integral_errors,err,error,*999)
      integral_error(1)=analytic_analysis_percentage_error(integral_errors(1,component_number),integral_errors(3,component_number))
      integral_error(2)=analytic_analysis_percentage_error(integral_errors(2,component_number),integral_errors(4,component_number))
      ghost_integral_error(1)=analytic_analysis_percentage_error(ghost_integral_errors(1,component_number), &
        & ghost_integral_errors(3,component_number))
      ghost_integral_error(2)=analytic_analysis_percentage_error(ghost_integral_errors(2,component_number), &
        & ghost_integral_errors(4,component_number))
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF 
    
    exits("AnalyticAnalysis_IntegralPercentageErrorGet")
    RETURN
999 errorsexits("AnalyticAnalysis_IntegralPercentageErrorGet",err,error)
    RETURN 1
    
  END SUBROUTINE analyticanalysis_integralpercentageerrorget
  
   !
  !================================================================================================================================
  !

  SUBROUTINE analyticanalysis_integralrelativeerrorget(FIELD,VARIABLE_TYPE,COMPONENT_NUMBER,INTEGRAL_ERROR, &
    & ghost_integral_error,err,error,*)
  
    !Argument variables   
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    REAL(DP), INTENT(OUT) :: INTEGRAL_ERROR(2)
    REAL(DP), INTENT(OUT) :: GHOST_INTEGRAL_ERROR(2)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
    REAL(DP), ALLOCATABLE :: INTEGRAL_ERRORS(:,:)
    REAL(DP), ALLOCATABLE :: GHOST_INTEGRAL_ERRORS(:,:)

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

    IF(ASSOCIATED(field)) THEN
      field_variable=>field%VARIABLE_TYPE_MAP(variable_type)%PTR
      CALL analytic_analysis_integral_errors(field_variable,integral_errors,ghost_integral_errors,err,error,*999)
      integral_error(1)=analytic_analysis_relative_error(integral_errors(1,component_number),integral_errors(3,component_number))
      integral_error(2)=analytic_analysis_relative_error(integral_errors(2,component_number),integral_errors(4,component_number))
      ghost_integral_error(1)=analytic_analysis_relative_error(ghost_integral_errors(1,component_number), &
        & ghost_integral_errors(3,component_number))
      ghost_integral_error(2)=analytic_analysis_relative_error(ghost_integral_errors(2,component_number), &
        & ghost_integral_errors(4,component_number))
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF 
    
    exits("AnalyticAnalysis_IntegralRelativeErrorGet")
    RETURN
999 errorsexits("AnalyticAnalysis_IntegralRelativeErrorGet",err,error)
    RETURN 1
    
  END SUBROUTINE analyticanalysis_integralrelativeerrorget

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

  SUBROUTINE analyticanalysis_absoluteerrorgetnode(FIELD,VARIABLE_TYPE,VERSION_NUMBER,DERIVATIVE_NUMBER,USER_NODE_NUMBER, &
    & component_number,VALUE,err,error,*)

    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: VERSION_NUMBER
    INTEGER(INTG), INTENT(IN) :: DERIVATIVE_NUMBER
    INTEGER(INTG), INTENT(IN) :: USER_NODE_NUMBER
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    REAL(DP), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    REAL(DP) :: NUMERICAL_VALUE, ANALYTIC_VALUE

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

    IF(ASSOCIATED(field)) THEN
      CALL field_parameter_set_get_node(field,variable_type,field_values_set_type,version_number,derivative_number, &
        & user_node_number,component_number,numerical_value,err,error,*999)
      CALL field_parameter_set_get_node(field,variable_type,field_analytic_values_set_type,version_number,derivative_number, &
        & user_node_number,component_number,analytic_value,err,error,*999)
      VALUE=analytic_analysis_absolute_error(numerical_value,analytic_value)
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF

    exits("AnalyticAnalysis_AbsoluteErrorGetNode")
    RETURN
999 errorsexits("AnalyticAnalysis_AbsoluteErrorGetNode",err,error)
    RETURN 1
  END SUBROUTINE analyticanalysis_absoluteerrorgetnode

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

  SUBROUTINE analyticanalysis_percentageerrorgetnode(FIELD,VARIABLE_TYPE,VERSION_NUMBER,DERIVATIVE_NUMBER,USER_NODE_NUMBER, &
    & component_number,VALUE,err,error,*)

    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: VERSION_NUMBER
    INTEGER(INTG), INTENT(IN) :: DERIVATIVE_NUMBER
    INTEGER(INTG), INTENT(IN) :: USER_NODE_NUMBER
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    REAL(DP), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    REAL(DP) :: NUMERICAL_VALUE, ANALYTIC_VALUE

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

    IF(ASSOCIATED(field)) THEN
      CALL field_parameter_set_get_node(field,variable_type,field_values_set_type,version_number,derivative_number, &
        & user_node_number,component_number,numerical_value,err,error,*999)
      CALL field_parameter_set_get_node(field,variable_type,field_analytic_values_set_type,version_number,derivative_number, &
        & user_node_number,component_number,analytic_value,err,error,*999)
      VALUE=analytic_analysis_percentage_error(numerical_value, analytic_value)
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF

    exits("AnalyticAnalysis_PercentageErrorGetNode")
    RETURN
999 errorsexits("AnalyticAnalysis_PercentageErrorGetNode",err,error)
    RETURN 1
  END SUBROUTINE analyticanalysis_percentageerrorgetnode



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

  SUBROUTINE analyticanalysis_relativeerrorgetnode(FIELD,VARIABLE_TYPE,VERSION_NUMBER,DERIVATIVE_NUMBER,USER_NODE_NUMBER, &
    & component_number,VALUE,err,error,*)

    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: VERSION_NUMBER
    INTEGER(INTG), INTENT(IN) :: DERIVATIVE_NUMBER
    INTEGER(INTG), INTENT(IN) :: USER_NODE_NUMBER
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    REAL(DP), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    REAL(DP) :: NUMERICAL_VALUE, ANALYTIC_VALUE

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

    IF(ASSOCIATED(field)) THEN
      CALL field_parameter_set_get_node(field,variable_type,field_values_set_type,version_number,derivative_number, &
        & user_node_number,component_number,numerical_value,err,error,*999)
      CALL field_parameter_set_get_node(field,variable_type,field_analytic_values_set_type,version_number,derivative_number, &
        & user_node_number,component_number,analytic_value,err,error,*999)
      VALUE=analytic_analysis_relative_error(numerical_value, analytic_value)
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF

    exits("AnalyticAnalysis_RelativeErrorGetNode")
    RETURN
999 errorsexits("AnalyticAnalysis_RelativeErrorGetNode",err,error)
    RETURN 1
  END SUBROUTINE analyticanalysis_relativeerrorgetnode

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

  SUBROUTINE analyticanalysis_absoluteerrorgetelement(FIELD,VARIABLE_TYPE,USER_ELEMENT_NUMBER,COMPONENT_NUMBER,VALUE,ERR, &
    & error,*)

    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: USER_ELEMENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    REAL(DP), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    REAL(DP) :: NUMERICAL_VALUE, ANALYTIC_VALUE

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

    IF(ASSOCIATED(field)) THEN
      CALL field_parameter_set_get_element(field,variable_type,field_values_set_type,user_element_number,component_number, &
        & numerical_value,err,error,*999)
      CALL field_parameter_set_get_element(field,variable_type,field_analytic_values_set_type,user_element_number, &
        & component_number,analytic_value,err,error,*999)
      VALUE=analytic_analysis_absolute_error(numerical_value,analytic_value)
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF

    exits("AnalyticAnalysis_AbsoluteErrorGetElement")
    RETURN
999 errorsexits("AnalyticAnalysis_AbsoluteErrorGetElement",err,error)
    RETURN 1
    
  END SUBROUTINE analyticanalysis_absoluteerrorgetelement

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

  SUBROUTINE analyticanalysis_percentageerrorgetelement(FIELD,VARIABLE_TYPE, &
    & user_element_number,component_number,VALUE,err,error,*)

    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: USER_ELEMENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    REAL(DP), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    REAL(DP) :: NUMERICAL_VALUE, ANALYTIC_VALUE

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

    IF(ASSOCIATED(field)) THEN
      CALL field_parameter_set_get_element(field,variable_type,field_values_set_type,user_element_number,component_number, &
        & numerical_value,err,error,*999)
      CALL field_parameter_set_get_element(field,variable_type,field_analytic_values_set_type,user_element_number, &
        & component_number,analytic_value,err,error,*999)
      VALUE=analytic_analysis_percentage_error(numerical_value, analytic_value)
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF

    exits("AnalyticAnalysis_PercentageErrorGetElement")
    RETURN
999 errorsexits("AnalyticAnalysis_PercentageErrorGetElement",err,error)
    RETURN 1
  END SUBROUTINE analyticanalysis_percentageerrorgetelement


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

  SUBROUTINE analyticanalysis_relativeerrorgetelement(FIELD,VARIABLE_TYPE,USER_ELEMENT_NUMBER,COMPONENT_NUMBER,VALUE,ERR, &
    & error,*)

    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: USER_ELEMENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    REAL(DP), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    REAL(DP) :: NUMERICAL_VALUE, ANALYTIC_VALUE

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

    IF(ASSOCIATED(field)) THEN
      CALL field_parameter_set_get_element(field,variable_type,field_values_set_type,user_element_number,component_number, &
        & numerical_value,err,error,*999)
      CALL field_parameter_set_get_element(field,variable_type,field_analytic_values_set_type,user_element_number, &
        & component_number,analytic_value,err,error,*999)
      VALUE=analytic_analysis_relative_error(numerical_value, analytic_value)
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF

    exits("AnalyticAnalysis_RelativeErrorGetElement")
    RETURN
999 errorsexits("AnalyticAnalysis_RelativeErrorGetElement",err,error)
    RETURN 1
    
  END SUBROUTINE analyticanalysis_relativeerrorgetelement

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

  SUBROUTINE analyticanalysis_absoluteerrorgetconstant(FIELD,VARIABLE_TYPE,COMPONENT_NUMBER,VALUE,ERR, &
    & error,*)

    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    REAL(DP), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    REAL(DP) :: NUMERICAL_VALUE, ANALYTIC_VALUE

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

    IF(ASSOCIATED(field)) THEN
      CALL field_parameter_set_get_constant(field,variable_type,field_values_set_type,component_number,numerical_value,err,error, &
        & *999)
      CALL field_parameter_set_get_constant(field,variable_type,field_analytic_values_set_type,component_number,analytic_value, &
        & err,error,*999)
      VALUE=analytic_analysis_absolute_error(numerical_value,analytic_value)
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF

    exits("AnalyticAnalysis_AbsoluteErrorGetConstant")
    RETURN
999 errorsexits("AnalyticAnalysis_AbsoluteErrorGetConstant",err,error)
    RETURN 1
    
  END SUBROUTINE analyticanalysis_absoluteerrorgetconstant

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

  SUBROUTINE analyticanalysis_percentageerrorgetconstant(FIELD,VARIABLE_TYPE,COMPONENT_NUMBER,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    REAL(DP), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    REAL(DP) :: NUMERICAL_VALUE, ANALYTIC_VALUE

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

    IF(ASSOCIATED(field)) THEN
      CALL field_parameter_set_get_constant(field,variable_type,field_values_set_type,component_number,numerical_value,err,error, &
        & *999)
      CALL field_parameter_set_get_constant(field,variable_type,field_analytic_values_set_type,component_number,analytic_value, &
        & err,error,*999)
      VALUE=analytic_analysis_percentage_error(numerical_value, analytic_value)
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF

    exits("AnalyticAnalysis_PercentageErrorGetConstant")
    RETURN
999 errorsexits("AnalyticAnalysis_PercentageErrorGetConstant",err,error)
    RETURN 1
    
  END SUBROUTINE analyticanalysis_percentageerrorgetconstant


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

  SUBROUTINE analyticanalysis_relativeerrorgetconstant(FIELD,VARIABLE_TYPE,COMPONENT_NUMBER,VALUE,ERR, &
    & error,*)

    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    REAL(DP), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    REAL(DP) :: NUMERICAL_VALUE, ANALYTIC_VALUE

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

    IF(ASSOCIATED(field)) THEN
      CALL field_parameter_set_get_constant(field,variable_type,field_values_set_type,component_number,numerical_value,err,error, &
        & *999)
      CALL field_parameter_set_get_constant(field,variable_type,field_analytic_values_set_type,component_number,analytic_value, &
        & err,error,*999)
      VALUE=analytic_analysis_relative_error(numerical_value, analytic_value)
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF

    exits("AnalyticAnalysis_RelativeErrorGetConstant")
    RETURN
999 errorsexits("AnalyticAnalysis_RelativeErrorGetConstant",err,error)
    RETURN 1
    
  END SUBROUTINE analyticanalysis_relativeerrorgetconstant

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

  SUBROUTINE analyticanalysis_rmserrorgetnode(FIELD,VARIABLE_TYPE,COMPONENT_NUMBER,ERROR_TYPE,LOCAL_RMS,LOCAL_GHOST_RMS, &
    & global_rms,err,error,*)
  
    !Argument variables   
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: ERROR_TYPE
    REAL(DP), INTENT(OUT) :: LOCAL_RMS(8)
    REAL(DP), INTENT(OUT) :: LOCAL_GHOST_RMS(8)
    REAL(DP), INTENT(OUT) :: GLOBAL_RMS(8)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    REAL(DP) :: ERROR_VALUE
    INTEGER(INTG) :: GHOST_NUMBER(8),NUMBER(8),MPI_IERROR
    REAL(DP) :: RMS_ERROR(8),GHOST_RMS_ERROR(8)
    TYPE(domain_nodes_type), POINTER :: NODES_DOMAIN
    INTEGER(INTG) :: node_idx,deriv_idx
        
    enters("ANALYTIC_ANALYSIS_RMS_PERCENTAGE_ERROR_GET_NODE",err,error,*999)

    IF(ASSOCIATED(field)) THEN
      nodes_domain=>field%VARIABLE_TYPE_MAP(variable_type)%PTR%COMPONENTS(component_number)%DOMAIN%TOPOLOGY%NODES
      IF(ASSOCIATED(nodes_domain)) THEN
        number=0
        rms_error=0.0_dp
        ghost_number=0
        ghost_rms_error=0.0_dp
        DO node_idx=1,nodes_domain%NUMBER_OF_NODES
          DO deriv_idx=1,nodes_domain%NODES(node_idx)%NUMBER_OF_DERIVATIVES
            SELECT CASE(error_type)
            CASE(absolute_error_type)
              !Default to version 1 of each node derivative
              CALL analyticanalysis_absoluteerrorgetnode(field,variable_type,1,deriv_idx,node_idx,component_number, &
                & error_value,err,error,*999)
            CASE(percentage_error_type)
              !Default to version 1 of each node derivative
              CALL analyticanalysis_percentageerrorgetnode(field,variable_type,1,deriv_idx,node_idx,component_number, &
                & error_value,err,error,*999)
            CASE(relative_error_type)
              !Default to version 1 of each node derivative
              CALL analyticanalysis_relativeerrorgetnode(field,variable_type,1,deriv_idx,node_idx,component_number, &
                & error_value,err,error,*999)
            CASE DEFAULT
              CALL flagerror("The error type is not valid!",err,error,*999)
            END SELECT
            !Accumlate the RMS errors
            number(deriv_idx)=number(deriv_idx)+1
            rms_error(deriv_idx)=rms_error(deriv_idx)+error_value*error_value
          ENDDO !deriv_idx
        ENDDO !node_idx
        DO node_idx=nodes_domain%NUMBER_OF_NODES+1,nodes_domain%TOTAL_NUMBER_OF_NODES
          DO deriv_idx=1,nodes_domain%NODES(node_idx)%NUMBER_OF_DERIVATIVES
            SELECT CASE(error_type)
            CASE(absolute_error_type)
              !Default to version 1 of each node derivative
              CALL analyticanalysis_absoluteerrorgetnode(field,variable_type,1,deriv_idx,node_idx,component_number, &
                & error_value,err,error,*999)
            CASE(percentage_error_type)
              !Default to version 1 of each node derivative
              CALL analyticanalysis_percentageerrorgetnode(field,variable_type,1,deriv_idx,node_idx,component_number, &
                & error_value,err,error,*999)
            CASE(relative_error_type)
              !Default to version 1 of each node derivative
              CALL analyticanalysis_relativeerrorgetnode(field,variable_type,1,deriv_idx,node_idx,component_number, &
                & error_value,err,error,*999)
            CASE DEFAULT
              CALL flagerror("The error type is not valid!",err,error,*999)
            END SELECT
            !Accumlate the RMS errors
            ghost_number(deriv_idx)=ghost_number(deriv_idx)+1
            ghost_rms_error(deriv_idx)=ghost_rms_error(deriv_idx)+error_value*error_value
          ENDDO !deriv_idx
        ENDDO !node_idx

        IF(computational_environment%NUMBER_COMPUTATIONAL_NODES>1) THEN
          IF(any(number>0)) THEN
            DO deriv_idx=1,8
              IF(number(deriv_idx)>0) THEN
                local_rms(deriv_idx)=sqrt(rms_error(deriv_idx)/number(deriv_idx))
              ENDIF
            ENDDO !deriv_idx
            DO deriv_idx=1,8
              IF(number(deriv_idx)>0) THEN
                local_ghost_rms(deriv_idx)=sqrt((rms_error(deriv_idx)+ghost_rms_error(deriv_idx))/(number(deriv_idx) &
                      & +ghost_number(deriv_idx)))
              ENDIF
            ENDDO !deriv_idx
            !Global RMS values
            !Collect the values across the ranks
            CALL mpi_allreduce(mpi_in_place,number,8,mpi_integer,mpi_sum,computational_environment%MPI_COMM,mpi_ierror)
            CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
            CALL mpi_allreduce(mpi_in_place,rms_error,8,mpi_double_precision,mpi_sum,computational_environment%MPI_COMM, &
              & mpi_ierror)
            CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
            DO deriv_idx=1,8
              IF(number(deriv_idx)>0) THEN
                global_rms(deriv_idx)=sqrt(rms_error(deriv_idx)/number(deriv_idx))
              ENDIF
            ENDDO !deriv_idx
          ENDIF
        ELSE
          IF(any(number>0)) THEN
            DO deriv_idx=1,8
              IF(number(deriv_idx)>0) THEN
                local_rms(deriv_idx)=sqrt(rms_error(deriv_idx)/number(deriv_idx))
                global_rms(deriv_idx)=local_rms(deriv_idx)
              ENDIF
            ENDDO !deriv_idx
          ENDIF
        ENDIF
      ELSE
        CALL flagerror("Nodes domain topology is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF 
    
    exits("AnalyticAnalysis_RMSErrorGetNode")
    RETURN
999 errorsexits("AnalyticAnalysis_RMSErrorGetNode",err,error)
    RETURN 1
  END SUBROUTINE analyticanalysis_rmserrorgetnode

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

  SUBROUTINE analyticanalysis_rmserrorgetelement(FIELD,VARIABLE_TYPE,COMPONENT_NUMBER,ERROR_TYPE,LOCAL_RMS,LOCAL_GHOST_RMS, &
    & global_rms,err,error,*)

    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(IN) :: VARIABLE_TYPE
    INTEGER(INTG), INTENT(IN) :: COMPONENT_NUMBER
    INTEGER(INTG), INTENT(IN) :: ERROR_TYPE
    REAL(DP), INTENT(OUT) :: LOCAL_RMS
    REAL(DP), INTENT(OUT) :: LOCAL_GHOST_RMS
    REAL(DP), INTENT(OUT) :: GLOBAL_RMS
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    REAL(DP) :: ERROR_VALUE
    INTEGER(INTG) :: GHOST_NUMBER,NUMBER,MPI_IERROR
    REAL(DP) :: RMS_ERROR,GHOST_RMS_ERROR
    TYPE(decomposition_type), POINTER :: DECOMPOSITION
    TYPE(decomposition_elements_type), POINTER :: ELEMENTS_DECOMPOSITION
    TYPE(decomposition_topology_type), POINTER :: DECOMPOSITION_TOPOLOGY
    TYPE(domain_type), POINTER :: DOMAIN
    TYPE(domain_elements_type), POINTER :: ELEMENTS_DOMAIN
    INTEGER(INTG) :: element_idx

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

    IF(ASSOCIATED(field)) THEN
      domain=>field%VARIABLE_TYPE_MAP(variable_type)%PTR%COMPONENTS(component_number)%DOMAIN
      elements_domain=>domain%TOPOLOGY%ELEMENTS
      IF(ASSOCIATED(elements_domain)) THEN
        decomposition=>domain%DECOMPOSITION
        IF(ASSOCIATED(decomposition)) THEN
          decomposition_topology=>decomposition%TOPOLOGY
          IF(ASSOCIATED(decomposition_topology)) THEN
            elements_decomposition=>decomposition_topology%ELEMENTS
            IF(ASSOCIATED(elements_decomposition)) THEN
              number=0
              rms_error=0.0_dp
              ghost_number=0
              ghost_rms_error=0.0_dp
              DO element_idx=1,elements_domain%NUMBER_OF_ELEMENTS
                SELECT CASE(error_type)
                CASE(absolute_error_type)
                  CALL analyticanalysis_absoluteerrorgetelement(field,variable_type,element_idx,component_number,error_value, &
                    & err,error,*999)
                CASE(percentage_error_type)
                  CALL analyticanalysis_percentageerrorgetelement(field,variable_type,element_idx,component_number, &
                    & error_value,err,error,*999)
                CASE(relative_error_type)
                  CALL analyticanalysis_relativeerrorgetelement(field,variable_type,element_idx,component_number,error_value, &
                    & err,error,*999)
                CASE DEFAULT
                  CALL flagerror("The error type is not valid!",err,error,*999)
                END SELECT
                number=number+1
                rms_error=rms_error+error_value*error_value
              ENDDO !element_idx
              DO element_idx=elements_domain%NUMBER_OF_ELEMENTS+1,elements_domain%TOTAL_NUMBER_OF_ELEMENTS
                SELECT CASE(error_type)
                CASE(absolute_error_type)
                  CALL analyticanalysis_absoluteerrorgetelement(field,variable_type,element_idx,component_number,error_value, &
                    & err,error,*999)
                CASE(percentage_error_type)
                  CALL analyticanalysis_percentageerrorgetelement(field,variable_type,element_idx,component_number, &
                    & error_value,err,error,*999)
                CASE(relative_error_type)
                  CALL analyticanalysis_relativeerrorgetelement(field,variable_type,element_idx,component_number,error_value, &
                    & err,error,*999)
                CASE DEFAULT
                  CALL flagerror("The error type is not valid!",err,error,*999)
                END SELECT
                ghost_number=ghost_number+1
                ghost_rms_error=ghost_rms_error+error_value*error_value
              ENDDO !element_idx
              IF(number>0) THEN
                IF(computational_environment%NUMBER_COMPUTATIONAL_NODES>1) THEN
                  !Local elements only
                  local_rms=sqrt(rms_error/number)
                  !Local and ghost elements
                  local_ghost_rms=sqrt((rms_error+ghost_rms_error)/(number+ghost_number))
                  !Global RMS values
                  !Collect the values across the ranks
                  CALL mpi_allreduce(mpi_in_place,number,1,mpi_integer,mpi_sum,computational_environment%MPI_COMM,mpi_ierror)
                  CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
                  CALL mpi_allreduce(mpi_in_place,rms_error,1,mpi_double_precision,mpi_sum,computational_environment%MPI_COMM, &
                    & mpi_ierror)
                  CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
                  global_rms=sqrt(rms_error/number)
                ENDIF
              ENDIF
            ELSE
              CALL flagerror("Decomposition topology elements is not associated.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Decomposition topology is not associated.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Domain decomposition is not associated.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Elements domain topology is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Field is not associated",err,error,*999)
    ENDIF

    exits("AnalyticAnalysis_RMSErrorGetElement")
    RETURN
999 errorsexits("AnalyticAnalysis_RMSErrorGetElement",err,error)
    RETURN 1
  END SUBROUTINE analyticanalysis_rmserrorgetelement

  !
  !================================================================================================================================
  !  
  
 
END MODULE analytic_analysis_routines