field_IO_routines.f90

Go to the documentation of this file.

MODULE field_io_routines
  USE base_routines
  USE lists
  USE basis_routines
  USE mesh_routines
  USE node_routines
  USE comp_environment
  USE coordinate_routines
  USE iso_varying_string
  USE machine_constants
  USE kinds
  USE field_routines
  USE iso_varying_string
  !USE, INTRINSIC :: ISO_C_BINDING
  USE strings
  USE types
  USE constants
#ifndef NOMPIMOD
  USE mpi
#endif
  USE cmiss_mpi
  USE input_output
  USE distributed_matrix_vector

#include "macros.h"  

  IMPLICIT NONE

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

  PRIVATE

  !Module parameters

  INTEGER(INTG), PARAMETER :: shape_size=3

  INTEGER(INTG), PARAMETER :: field_io_field_label=1
  INTEGER(INTG), PARAMETER :: field_io_variable_label=2
  INTEGER(INTG), PARAMETER :: field_io_component_label=3
  INTEGER(INTG), PARAMETER :: field_io_derivative_label=4

  INTEGER(INTG), PARAMETER :: field_io_scale_factors_number_type=5
  INTEGER(INTG), PARAMETER :: field_io_scale_factors_property_type=6

  !Module types

  TYPE mesh_elements_type_ptr_type
    TYPE(meshelementstype), POINTER :: ptr
  END TYPE mesh_elements_type_ptr_type

  TYPE field_variable_component_ptr_type
    TYPE(field_variable_component_type), POINTER :: ptr
  END TYPE field_variable_component_ptr_type

  TYPE field_io_component_info_set
    LOGICAL :: same_header
    INTEGER(INTG) :: number_of_components
    !attention: the pointers in COMPONENTS(:) point to those nodal components which are in the same local domain in current implementation
    !it may be replaced in the future implementation
    TYPE(field_variable_component_ptr_type), ALLOCATABLE:: components(:)
    INTEGER(INTG), ALLOCATABLE:: component_versions(:)
  END TYPE field_io_component_info_set

  TYPE field_io_component_info_set_ptr_type
    TYPE(field_io_component_info_set), POINTER :: ptr
  END TYPE field_io_component_info_set_ptr_type

  TYPE field_io_info_set
    TYPE(fields_type), POINTER :: fields
    INTEGER(INTG) :: number_of_entries
    !Interesting thing: pointer here, also means dymanically allocated attibute
    INTEGER(INTG), ALLOCATABLE:: list_of_global_number(:)
    TYPE(field_io_component_info_set_ptr_type), ALLOCATABLE:: component_info_set(:)
  END TYPE field_io_info_set

  !Module variables

  !Interfaces
  INTERFACE
    FUNCTION fieldexport_opensession( exportType, filename, handle ) &
      & bind(c,name="FieldExport_OpenSession")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: exportType
      CHARACTER(C_CHAR), INTENT(IN) :: filename(*)
      INTEGER(C_INT), INTENT(OUT) :: handle
      INTEGER(C_INT) :: FieldExport_OpenSession
    END FUNCTION fieldexport_opensession

    FUNCTION fieldexport_group( handle, groupName ) &
      & bind(c,name="FieldExport_Group")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      CHARACTER(C_CHAR), INTENT(IN) :: groupName(*)
      INTEGER(C_INT) :: FieldExport_Group
    END FUNCTION fieldexport_group

    FUNCTION fieldexport_meshdimensions( handle, meshDimensions , basisType) &
      & bind(c,name="FieldExport_MeshDimensions")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: meshDimensions
      INTEGER(C_INT), VALUE :: basisType
      INTEGER(C_INT) :: FieldExport_MeshDimensions
    END FUNCTION fieldexport_meshdimensions

    FUNCTION fieldexport_scalingfactorcount( handle, scalingFactorCount ) &
      & bind(c,name="FieldExport_ScalingFactorCount")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: scalingFactorCount
      INTEGER(C_INT) :: FieldExport_ScalingFactorCount
    END FUNCTION fieldexport_scalingfactorcount

    FUNCTION fieldexport_scalefactors( handle, numberOfXi, interpolationXi, numberOfScaleFactors ) &
      & bind(c,name="FieldExport_ScaleFactors")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: numberOfXi
      TYPE(c_ptr), VALUE :: interpolationXi
      INTEGER(C_INT), VALUE :: numberOfScaleFactors
      INTEGER(C_INT) :: FieldExport_ScaleFactors
    END FUNCTION fieldexport_scalefactors

    FUNCTION fieldexport_nodecount( handle, nodeCount ) &
      & bind(c,name="FieldExport_NodeCount")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: nodeCount
      INTEGER(C_INT) :: FieldExport_NodeCount
    END FUNCTION fieldexport_nodecount

    FUNCTION fieldexport_fieldcount( handle, fieldCount ) &
      & bind(c,name="FieldExport_FieldCount")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: fieldCount
      INTEGER(C_INT) :: FieldExport_FieldCount
    END FUNCTION fieldexport_fieldcount

    FUNCTION fieldexport_coordinatevariable( handle, variableName, variableNumber, coordinateSystemType, componentCount ) &
      & bind(c,name="FieldExport_CoordinateVariable")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      CHARACTER(LEN=1, KIND=C_CHAR) :: variableName(*)
      INTEGER(C_INT), VALUE :: variableNumber
      INTEGER(C_INT), VALUE :: coordinateSystemType
      INTEGER(C_INT), VALUE :: componentCount
      INTEGER(C_INT) :: FieldExport_CoordinateVariable
    END FUNCTION fieldexport_coordinatevariable

    FUNCTION fieldexport_variable( handle, variableName, variableNumber, fieldType, variableType, componentCount ) &
      & bind(c,name="FieldExport_Variable")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      CHARACTER(LEN=1, KIND=C_CHAR) :: variableName(*)
      INTEGER(C_INT), VALUE :: variableNumber
      INTEGER(C_INT), VALUE :: fieldType
      INTEGER(C_INT), VALUE :: variableType
      INTEGER(C_INT), VALUE :: componentCount
      INTEGER(C_INT) :: FieldExport_Variable
    END FUNCTION fieldexport_variable

    FUNCTION fieldexport_coordinatecomponent( handle, coordinateSystemType, componentNumber, interpType, &
      & numberofxi, interpolationxi ) &
      & bind(c,name="FieldExport_CoordinateComponent")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: coordinateSystemType
      INTEGER(C_INT), VALUE :: componentNumber
      INTEGER(C_INT), VALUE :: interpType
      INTEGER(C_INT), VALUE :: numberOfXi
      TYPE(c_ptr), VALUE :: interpolationXi
      INTEGER(C_INT) :: FieldExport_CoordinateComponent
    END FUNCTION fieldexport_coordinatecomponent

    FUNCTION fieldexport_component( handle, componentNumber, interpType, numberOfXi, interpolationXi ) &
      & bind(c,name="FieldExport_Component")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: componentNumber
      INTEGER(C_INT), VALUE :: interpType
      INTEGER(C_INT), VALUE :: numberOfXi
      TYPE(c_ptr), VALUE :: interpolationXi
      INTEGER(C_INT) :: FieldExport_Component
    END FUNCTION fieldexport_component

    FUNCTION fieldexport_elementgridsize( handle, headerType, numberOfXi, numberGauss ) &
      & bind(c,name="FieldExport_ElementGridSize")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: headerType
      INTEGER(C_INT), VALUE :: numberOfXi
      TYPE(c_ptr), VALUE :: numberGauss
      INTEGER(C_INT) :: FieldExport_ElementGridSize
    END FUNCTION fieldexport_elementgridsize

    FUNCTION fieldexport_nodescaleindexes( handle, nodeCount, derivativeCount, elementDerivatives, nodeIndexes, &
      & scaleindexes ) &
      & bind(c,name="FieldExport_NodeScaleIndexes")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: nodeCount
      TYPE(c_ptr), VALUE :: derivativeCount
      TYPE(c_ptr), VALUE :: elementDerivatives
      TYPE(c_ptr), VALUE :: nodeIndexes
      TYPE(c_ptr), VALUE :: scaleIndexes
      INTEGER(C_INT) :: FieldExport_NodeScaleIndexes
    END FUNCTION fieldexport_nodescaleindexes

    FUNCTION fieldexport_elementindex( handle, dimensionCount, elementIndex ) &
      & bind(c,name="FieldExport_ElementIndex")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: dimensionCount
      INTEGER(C_INT), VALUE :: elementIndex
      INTEGER(C_INT) :: FieldExport_ElementIndex
    END FUNCTION fieldexport_elementindex

    FUNCTION fieldexport_elementnodeindices( handle, nodeCount, nodeIndices ) &
      & bind(c,name="FieldExport_ElementNodeIndices")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: nodeCount
      TYPE(c_ptr), VALUE :: nodeIndices
      INTEGER(C_INT) :: FieldExport_ElementNodeIndices
    END FUNCTION fieldexport_elementnodeindices

    FUNCTION fieldexport_elementnodescales( handle, isFirstSet, scaleCount, scales ) &
      & bind(c,name="FieldExport_ElementNodeScales")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: isFirstSet
      INTEGER(C_INT), VALUE :: scaleCount
      TYPE(c_ptr), VALUE :: scales
      INTEGER(C_INT) :: FieldExport_ElementNodeScales
    END FUNCTION fieldexport_elementnodescales

    FUNCTION fieldexport_elementgridvalues( handle, isFirstSet, numberOfXi, elementValue ) &
      & bind(c,name="FieldExport_ElementGridValues")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: isFirstSet
      INTEGER(C_INT), VALUE :: numberOfXi
      REAL(DP), VALUE :: elementValue
      INTEGER(C_INT) :: FieldExport_ElementGridValues
    END FUNCTION fieldexport_elementgridvalues

    FUNCTION fieldexport_nodevalues( handle, nodeNumber, valueCount, nodeValues ) &
      & bind(c,name="FieldExport_NodeValues")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: nodeNumber
      INTEGER(C_INT), VALUE :: valueCount
      TYPE(c_ptr), VALUE :: nodeValues
      INTEGER(C_INT) :: FieldExport_NodeValues
    END FUNCTION fieldexport_nodevalues

    FUNCTION fieldexport_closesession( handle ) &
      & bind(c,name="FieldExport_CloseSession")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT) :: FieldExport_CloseSession
    END FUNCTION fieldexport_closesession

    FUNCTION fieldexport_coordinatederivativeindices( handle, componentNumber, coordinateSystemType, numberOfDerivatives,  &
      & derivatives, valueindex ) bind(c,name="FieldExport_CoordinateDerivativeIndices")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: componentNumber
      INTEGER(C_INT), VALUE :: coordinateSystemType
      INTEGER(C_INT), VALUE :: numberOfDerivatives
      TYPE(c_ptr), VALUE :: derivatives
      INTEGER(C_INT), VALUE :: valueIndex
      INTEGER(C_INT) :: FieldExport_CoordinateDerivativeIndices
    END FUNCTION fieldexport_coordinatederivativeindices

    FUNCTION fieldexport_derivativeindices( handle, componentNumber, fieldType, variableType, numberOfDerivatives, &
      & derivatives, valueindex ) bind(c,name="FieldExport_DerivativeIndices")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: componentNumber
      INTEGER(C_INT), VALUE :: fieldType
      INTEGER(C_INT), VALUE :: variableType
      INTEGER(C_INT), VALUE :: numberOfDerivatives
      TYPE(c_ptr), VALUE :: derivatives
      INTEGER(C_INT), VALUE :: valueIndex
      INTEGER(C_INT) :: FieldExport_DerivativeIndices
    END FUNCTION fieldexport_derivativeindices

    FUNCTION fieldexport_endcomponent(handle) BIND(C,NAME="FieldExport_EndComponent")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT) :: FieldExport_EndComponent
    END FUNCTION fieldexport_endcomponent

    FUNCTION fieldexport_versioninfo(handle, numberOfVersions) BIND(C,NAME="FieldExport_VersionInfo")
      USE types
      USE iso_c_binding
      INTEGER(C_INT), VALUE :: handle
      INTEGER(C_INT), VALUE :: numberOfVersions
      INTEGER(C_INT) :: FieldExport_VersionInfo
    END FUNCTION fieldexport_versioninfo

  END INTERFACE

  INTERFACE reallocate
    MODULE PROCEDURE reallocate_int
    MODULE PROCEDURE reallocate_real
    MODULE PROCEDURE reallocate_string
    MODULE PROCEDURE reallocate_elements
    MODULE PROCEDURE reallocate_components
    MODULE PROCEDURE reallocate_basis
    MODULE PROCEDURE reallocate_field
  END INTERFACE !REALLOCATE

  INTERFACE grow_array
    MODULE PROCEDURE grow_array_int
    MODULE PROCEDURE grow_array_real
    MODULE PROCEDURE grow_array_components
  END INTERFACE !GROW_ARRAY

  INTERFACE checked_deallocate
    MODULE PROCEDURE checked_deallocate_int
    MODULE PROCEDURE checked_deallocate_real
    MODULE PROCEDURE checked_deallocate_2d_int
    MODULE PROCEDURE checked_deallocate_str
    MODULE PROCEDURE checked_deallocate_components
    MODULE PROCEDURE checked_deallocate_elements
    MODULE PROCEDURE checked_deallocate_basis
    MODULE PROCEDURE checked_deallocate_field
  END INTERFACE !CHECKED_DEALLOCATE

  PUBLIC :: field_io_fields_import, field_io_nodes_export, field_io_elements_export


CONTAINS

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

  SUBROUTINE reallocate_int( array, newSize, errorMessage, ERR, ERROR, * )
    INTEGER(INTG), ALLOCATABLE, INTENT(INOUT) :: array(:)
    INTEGER(INTG), INTENT(IN) :: newSize
    CHARACTER(LEN=*), INTENT(IN) :: errorMessage
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

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

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

    ALLOCATE( array( newsize ), stat = err )
    IF( err /= 0 ) CALL flagerror( errormessage, err, error, *999)

    array(:) = 0

    exits("REALLOCATE_INT")
    RETURN
999 errorsexits("REALLOCATE_INT",err,error)
    RETURN 1
  END SUBROUTINE reallocate_int

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

  SUBROUTINE reallocate_real( array, newSize, errorMessage, ERR, ERROR, * )
    REAL(DP), ALLOCATABLE, INTENT(INOUT) :: array(:)
    INTEGER(INTG), INTENT(IN) :: newSize
    CHARACTER(LEN=*), INTENT(IN) :: errorMessage
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

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

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

    ALLOCATE( array( newsize ), stat = err )
    IF( err /= 0 ) CALL flagerror( errormessage, err, error, *999)

    array(:) = 0

    exits("REALLOCATE_REAL")
    RETURN
999 errorsexits("REALLOCATE_REAL",err,error)
    RETURN 1
  END SUBROUTINE reallocate_real

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

  SUBROUTINE reallocate_string( array, newSize, errorMessage, ERR, ERROR, * )
    TYPE(varying_string), ALLOCATABLE, INTENT(INOUT) :: array(:)
    INTEGER(INTG), INTENT(IN) :: newSize
    CHARACTER(LEN=*), INTENT(IN) :: errorMessage
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

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

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

    ALLOCATE( array( newsize ), stat = err )
    IF( err /= 0 ) CALL flagerror( errormessage, err, error, *999)

    exits("REALLOCATE_STRING")
    RETURN
999 errorsexits("REALLOCATE_STRING",err,error)
    RETURN 1
  END SUBROUTINE reallocate_string

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

  SUBROUTINE reallocate_components( array, newSize, errorMessage, ERR, ERROR, * )
    TYPE(field_variable_component_ptr_type), ALLOCATABLE, INTENT(INOUT) :: array(:)
    INTEGER(INTG), INTENT(IN) :: newSize
    CHARACTER(LEN=*), INTENT(IN) :: errorMessage
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

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

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

    ALLOCATE( array( newsize ), stat = err )
    IF( err /= 0 ) CALL flagerror( errormessage, err, error, *999)

    exits("REALLOCATE_COMPONENTS")
    RETURN
999 errorsexits("REALLOCATE_COMPONENTS",err,error)
    RETURN 1
  END SUBROUTINE reallocate_components

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

  SUBROUTINE reallocate_basis( array, newSize, errorMessage, ERR, ERROR, * )
    TYPE(basis_ptr_type), ALLOCATABLE, INTENT(INOUT) :: array(:)
    INTEGER(INTG), INTENT(IN) :: newSize
    CHARACTER(LEN=*), INTENT(IN) :: errorMessage
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

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

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

    ALLOCATE( array( newsize ), stat = err )
    IF( err /= 0 ) CALL flagerror( errormessage, err, error, *999)

    exits("REALLOCATE_BASIS")
    RETURN
999 errorsexits("REALLOCATE_BASIS",err,error)
    RETURN 1
  END SUBROUTINE reallocate_basis

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

  SUBROUTINE reallocate_field( array, newSize, errorMessage, ERR, ERROR, * )
    TYPE(field_ptr_type), ALLOCATABLE, INTENT(INOUT) :: array(:)
    INTEGER(INTG), INTENT(IN) :: newSize
    CHARACTER(LEN=*), INTENT(IN) :: errorMessage
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

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

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

    ALLOCATE( array( newsize ), stat = err )
    IF( err /= 0 ) CALL flagerror( errormessage, err, error, *999)

    exits("REALLOCATE_FIELD")
    RETURN
999 errorsexits("REALLOCATE_FIELD",err,error)
    RETURN 1
  END SUBROUTINE reallocate_field

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

  SUBROUTINE reallocate_elements( array, newSize, errorMessage, ERR, ERROR, * )
    TYPE(mesh_elements_type_ptr_type), ALLOCATABLE, INTENT(INOUT) :: array(:)
    INTEGER(INTG), INTENT(IN) :: newSize
    CHARACTER(LEN=*), INTENT(IN) :: errorMessage
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

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

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

    ALLOCATE( array( newsize ), stat = err )
    IF( err /= 0 ) CALL flagerror( errormessage, err, error, *999)

    exits("REALLOCATE_ELEMENTS")
    RETURN
999 errorsexits("REALLOCATE_ELEMENTS",err,error)
    RETURN 1
  END SUBROUTINE reallocate_elements

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

  SUBROUTINE reallocate_2d( array, newSize1, newSize2, errorMessage, ERR, ERROR, * )
    INTEGER(INTG), ALLOCATABLE, INTENT(INOUT) :: array(:,:)
    INTEGER(INTG), INTENT(IN) :: newSize1
    INTEGER(INTG), INTENT(IN) :: newSize2
    CHARACTER(LEN=*), INTENT(IN) :: errorMessage
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

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

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

    ALLOCATE( array( newsize1, newsize2 ), stat = err )
    IF( err /= 0 ) CALL flagerror( errormessage, err, error, *999)

    array(:,:) = 0

    exits("REALLOCATE_2D")
    RETURN
999 errorsexits("REALLOCATE_2D",err,error)
    RETURN 1
  END SUBROUTINE reallocate_2d

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

  SUBROUTINE grow_array_int( array, delta, errorMessage, ERR, ERROR, * )
    INTEGER(INTG), ALLOCATABLE, INTENT(INOUT) :: array(:)
    INTEGER(INTG), INTENT(IN) :: delta
    CHARACTER(LEN=*), INTENT(IN) :: errorMessage
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    INTEGER(INTG), ALLOCATABLE :: tempArray(:)
    INTEGER(INTG) :: oldSize

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

    IF( .NOT.ALLOCATED( array ) ) THEN
      CALL reallocate( array, delta, errormessage, err, error, *999 )
      RETURN
    ENDIF

    oldsize = SIZE( array )

    CALL reallocate( temparray, oldsize, errormessage, err, error, *999 )

    temparray(:) = array(:)

    CALL reallocate( array, oldsize + delta, errormessage, err, error, *999 )

    array(1:oldsize) = temparray(:)

    DEALLOCATE( temparray )

    exits("GROW_ARRAY_INT")
    RETURN
999 errorsexits("GROW_ARRAY_INT",err,error)
    RETURN 1
  END SUBROUTINE grow_array_int

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

  SUBROUTINE grow_array_real( array, delta, errorMessage, ERR, ERROR, * )
    REAL(C_DOUBLE), ALLOCATABLE, INTENT(INOUT) :: array(:)
    INTEGER(INTG), INTENT(IN) :: delta
    CHARACTER(LEN=*), INTENT(IN) :: errorMessage
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    REAL(C_DOUBLE), ALLOCATABLE :: tempArray(:)
    INTEGER(INTG) :: oldSize

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

    IF( .NOT.ALLOCATED( array ) ) THEN
      CALL reallocate( array, delta, errormessage, err, error, *999 )
      RETURN
    ENDIF

    oldsize = SIZE( array )

    CALL reallocate( temparray, oldsize, errormessage, err, error, *999 )

    temparray(:) = array(:)

    CALL reallocate( array, oldsize + delta, errormessage, err, error, *999 )

    array(1:oldsize) = temparray(:)

    DEALLOCATE( temparray )

    exits("GROW_ARRAY_REAL")
    RETURN
999 errorsexits("GROW_ARRAY_REAL",err,error)
    RETURN 1
  END SUBROUTINE grow_array_real

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

  SUBROUTINE grow_array_components( array, delta, errorMessage, ERR, ERROR, * )
    TYPE(field_variable_component_ptr_type), ALLOCATABLE, INTENT(INOUT) :: array(:)
    INTEGER(INTG), INTENT(IN) :: delta
    CHARACTER(LEN=*), INTENT(IN) :: errorMessage
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    TYPE(field_variable_component_ptr_type), ALLOCATABLE :: tempArray(:)
    INTEGER(INTG) :: oldSize

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

    IF( .NOT.ALLOCATED( array ) ) THEN
      CALL reallocate( array, delta, errormessage, err, error, *999 )
      RETURN
    ENDIF

    oldsize = SIZE( array )

    CALL reallocate( temparray, oldsize, errormessage, err, error, *999 )

    temparray(:) = array(:)

    CALL reallocate( array, oldsize + delta, errormessage, err, error, *999 )

    array(1:oldsize) = temparray(:)

    DEALLOCATE( temparray )

    exits("GROW_ARRAY_COMPONENTS")
    RETURN
999 errorsexits("GROW_ARRAY_COMPONENTS",err,error)
    RETURN 1
  END SUBROUTINE grow_array_components

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

  SUBROUTINE checked_deallocate_int( array )
    INTEGER(INTG), ALLOCATABLE, INTENT(INOUT) :: array(:)

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

  END SUBROUTINE checked_deallocate_int

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

  SUBROUTINE checked_deallocate_real( array )
    REAL(DP), ALLOCATABLE, INTENT(INOUT) :: array(:)

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

  END SUBROUTINE checked_deallocate_real

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

  SUBROUTINE checked_deallocate_2d_int( array )
    INTEGER(INTG), ALLOCATABLE, INTENT(INOUT) :: array(:,:)

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

  END SUBROUTINE checked_deallocate_2d_int

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

  SUBROUTINE checked_deallocate_components( array )
    TYPE(field_variable_component_ptr_type), ALLOCATABLE, INTENT(INOUT) :: array(:)

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

  END SUBROUTINE checked_deallocate_components

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

  SUBROUTINE checked_deallocate_str( array )
    TYPE(varying_string), ALLOCATABLE, INTENT(INOUT) :: array(:)

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

  END SUBROUTINE checked_deallocate_str

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

  SUBROUTINE checked_deallocate_elements( array )
    TYPE(mesh_elements_type_ptr_type), ALLOCATABLE, INTENT(INOUT) :: array(:)

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

  END SUBROUTINE checked_deallocate_elements

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

  SUBROUTINE checked_deallocate_basis( array )
    TYPE(basis_ptr_type), ALLOCATABLE, INTENT(INOUT) :: array(:)

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

  END SUBROUTINE checked_deallocate_basis

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

  SUBROUTINE checked_deallocate_field( array )
    TYPE(field_ptr_type), ALLOCATABLE, INTENT(INOUT) :: array(:)

    IF( ALLOCATED( array ) ) THEN
      DEALLOCATE( array )
    ENDIF

  END SUBROUTINE checked_deallocate_field

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

  SUBROUTINE field_io_field_info(STRING, LABEL_TYPE, FIELD_TYPE , ERR, ERROR, *)
    !Argument variables
    TYPE(varying_string), INTENT(IN) :: STRING
    INTEGER(INTG), INTENT(IN) :: LABEL_TYPE
    INTEGER(INTG), INTENT(INOUT) :: FIELD_TYPE
    !REAL(DP), OPTIONAL, INTENT(INOUT) :: FOCUS
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: pos
    TYPE(varying_string) :: LINE, KEYWORD

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

    line=string

    SELECT CASE(label_type)
      CASE(field_io_field_label)
        pos=index(line, ",")
        line=remove(line, 1, pos)
        pos=index(line, ",")
        keyword=extract(line, 1, pos-1)
        line=remove(line, 1,pos)
        keyword=adjustl(keyword)
        keyword=trim(keyword)
        IF(keyword=="coordinate") THEN
          field_type=field_geometric_type
        ELSE IF (keyword=="anatomical") THEN
          field_type=field_fibre_type
        ELSE
          field_type=-1
          CALL flagerror("Cannot find corresponding field type from input string",err,error,*999)
        ENDIF
      CASE DEFAULT
        CALL flagerror("Cannot find any information from input string",err,error,*999)
    END SELECT !CASE(LABEL_TYPE)

    exits("FIELD_IO_FIELD_INFO")
    RETURN
999 errorsexits("FIELD_IO_FIELD_INFO",err,error)
    RETURN 1
  END SUBROUTINE field_io_field_info


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

  FUNCTION field_io_derivative_info(LINE, ERR, ERROR)
    !Argument variables
    TYPE(varying_string), INTENT(IN) :: LINE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) ::FIELD_IO_DERIVATIVE_INFO

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

    IF("d/ds1"==line) THEN
      field_io_derivative_info=part_deriv_s1
    ELSE IF("d2/ds1ds1"==line) THEN
      field_io_derivative_info=part_deriv_s1_s1
    ELSE IF("d/ds2"==line) THEN
      field_io_derivative_info=part_deriv_s2
    ELSE IF("d2/ds2ds2"==line) THEN
      field_io_derivative_info=part_deriv_s2_s2
    ELSE IF("d/ds3"==line) THEN
      field_io_derivative_info=part_deriv_s1_s2
    ELSE IF("d2/ds3ds3"==line) THEN
      field_io_derivative_info=part_deriv_s3
    ELSE IF("d2/ds3ds3"==line) THEN
      field_io_derivative_info=part_deriv_s3_s3
    ELSE IF("d2/ds1ds3"==line) THEN
      field_io_derivative_info=part_deriv_s1_s3
    ELSE IF("d2/ds2ds3"==line) THEN
      field_io_derivative_info=part_deriv_s2_s3
    ELSE IF("d3/ds1ds2ds3"==line) THEN
      field_io_derivative_info=part_deriv_s1_s2_s3
    ELSE IF("d/ds4"==line) THEN
      field_io_derivative_info=part_deriv_s4
    ELSE IF("d2/ds4ds4"==line) THEN
      field_io_derivative_info=part_deriv_s4_s4
    ELSE IF("d2/ds1ds4"==line) THEN
      field_io_derivative_info=part_deriv_s1_s4
    ELSE IF("d2/ds2ds4"==line) THEN
      field_io_derivative_info=part_deriv_s2_s4
    ELSE IF("d2/ds3ds4"==line) THEN
      field_io_derivative_info=part_deriv_s3_s4
    ELSE IF("d3/ds1ds2ds4"==line) THEN
      field_io_derivative_info=part_deriv_s1_s2_s4
    ELSE IF("d3/ds1ds3ds4"==line) THEN
      field_io_derivative_info=part_deriv_s1_s3_s4
    ELSE IF("d3/ds2ds3ds4"==line) THEN
      field_io_derivative_info=part_deriv_s2_s3_s4
    ELSE IF("d3/ds1ds4ds4"==line) THEN
      field_io_derivative_info=part_deriv_s1_s4_s4
    ELSE IF("d3/ds2ds4ds4"==line) THEN
      field_io_derivative_info=part_deriv_s2_s4_s4
    ELSE IF("d3/ds3ds4ds4"==line) THEN
      field_io_derivative_info=part_deriv_s3_s4_s4
    ELSE IF("d3/ds4ds4ds4"==line) THEN
      field_io_derivative_info=part_deriv_s4_s4_s4
    ELSE
      field_io_derivative_info=-1
      CALL flagerror("Could not recognize derivatives from input string",err,error,*999)
    ENDIF

    exits("FIELD_IO_DERIVATIVE_INFO")
    RETURN
999 errorsexits("FIELD_IO_DERIVATIVE_INFO",err,error)
  END FUNCTION field_io_derivative_info

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

  FUNCTION field_io_element_derivative_index(ELEMENT, DERIVATIVE_NUMBER, NODE_NUMBER, ERR, ERROR)
    !Argument variables
    TYPE(domain_element_type), INTENT(IN) :: ELEMENT
    INTEGER(INTG), INTENT(IN) :: DERIVATIVE_NUMBER
    INTEGER(INTG), INTENT(IN) :: NODE_NUMBER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Return variable
    INTEGER(INTG) :: FIELD_IO_ELEMENT_DERIVATIVE_INDEX
    !Local Variables
    INTEGER(INTG) :: VERSION_NUMBER,NUMBER_OF_DERIVATIVES

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

    version_number=element%elementVersions(derivative_number, node_number)
    number_of_derivatives=element%BASIS%NUMBER_OF_DERIVATIVES(node_number)
    field_io_element_derivative_index=(version_number-1)*number_of_derivatives + &
      & element%ELEMENT_DERIVATIVES(derivative_number, node_number)

    exits("FIELD_IO_ELEMENT_DERIVATIVE_INDEX")
    RETURN
999 errorsexits("FIELD_IO_ELEMENT_DERIVATIVE_INDEX",err,error)
  END FUNCTION field_io_element_derivative_index

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

  SUBROUTINE field_io_create_fields(NAME, REGION, DECOMPOSITION, FIELD_VALUES_SET_TYPE, NUMBER_OF_FIELDS, &
    !&USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER,
    &mesh_components_of_field_components, components_in_fields, number_of_exnode_files, &
    &master_computational_number, my_computational_node_number, field_scaling_type, err, error, *)
    !Argument variables
    TYPE(varying_string), INTENT(IN) :: NAME
    TYPE(region_type), POINTER :: REGION
    TYPE(decomposition_type), POINTER :: DECOMPOSITION
    INTEGER(INTG), INTENT(IN) :: FIELD_VALUES_SET_TYPE
    INTEGER(INTG), INTENT(IN) :: NUMBER_OF_FIELDS
    !INTEGER(INTG), INTENT(IN) :: USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER(:)
    INTEGER(INTG), INTENT(IN) :: MESH_COMPONENTS_OF_FIELD_COMPONENTS(:)
    INTEGER(INTG), INTENT(IN) :: COMPONENTS_IN_FIELDS(:)
    INTEGER(INTG), INTENT(IN) :: NUMBER_OF_EXNODE_FILES
    INTEGER(INTG), INTENT(IN) :: MASTER_COMPUTATIONAL_NUMBER
    INTEGER(INTG), INTENT(IN) :: my_computational_node_number
    INTEGER(INTG), INTENT(IN) :: FIELD_SCALING_TYPE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(field_type), POINTER :: FIELD
    TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES
    TYPE(varying_string), ALLOCATABLE :: LIST_STR(:)
    TYPE(varying_string) :: FILE_NAME, FILE_STATUS, LINE, LINE1
    TYPE(varying_string) :: CMISS_KEYWORD_FIELDS, CMISS_KEYWORD_NODE, CMISS_KEYWORD_COMPONENTS
    TYPE(varying_string) :: CMISS_KEYWORD_VALUE_INDEX, CMISS_KEYWORD_DERIVATIVE
    INTEGER(INTG), ALLOCATABLE :: tmp_pointer(:), LIST_DEV(:), LIST_DEV_POS(:)
    INTEGER(INTG) :: FILE_ID
    !INTEGER(INTG) :: NUMBER_FIELDS
    INTEGER(INTG) :: NODAL_USER_NUMBER, NODAL_LOCAL_NUMBER, FIELDTYPE, NUMBER_NODAL_VALUE_LINES, NUMBER_OF_LINES, &
      & NUMBER_OF_COMPONENTS !, LABEL_TYPE, FOCUS
    INTEGER(INTG) :: MPI_IERROR
    INTEGER(INTG) :: idx_comp, idx_comp1, pos, idx_field, idx_exnode, idx_nodal_line, idx_node
    INTEGER(INTG) :: idx_variable, idx_dev, idx_dev1, total_number_of_comps, total_number_of_devs, number_of_devs !idx_variable1
    INTEGER(INTG) :: number_of_comps, VARIABLE_IDX,variable_type
    REAL(DP), ALLOCATABLE :: LIST_DEV_VALUE(:)
    LOGICAL :: SECTION_START, FILE_END, NODE_SECTION, FILE_OPEN, NODE_IN_DOMAIN


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

    IF(.NOT.ASSOCIATED(decomposition)) THEN
      CALL flagerror("decomposition is NOT associated before importing data",err,error,*999)
      GOTO 999
    ENDIF

    IF(.NOT.ASSOCIATED(region)) THEN
      CALL flagerror("region is NOT associated before importing data",err,error,*999)
      GOTO 999
    ENDIF

    cmiss_keyword_fields="#Fields="
    cmiss_keyword_components="#Components="
    cmiss_keyword_value_index="Value index="
    cmiss_keyword_derivative="#Derivatives="
    cmiss_keyword_node="Node:"

    file_end=.false.
    idx_exnode=0
    file_id=1030
    file_status="OLD"
    total_number_of_comps=0
    number_nodal_value_lines=5
    number_of_components=sum(components_in_fields)

    !checking the field strings in exnode files
    IF(master_computational_number==my_computational_node_number) THEN

      CALL reallocate( list_str, number_of_fields, "can not allocate list of strings for fields", err, error, *999 )

      DO WHILE(idx_exnode<number_of_exnode_files)

        file_id=1030+idx_exnode
        !checking the next file
        file_name=name//".part"//trim(number_to_vstring(idx_exnode,"*",err,error))//".exnode"
        !INQUIRE(FILE=CHAR(FILE_NAME), OPENED=FILE_OPEN)
        CALL field_io_fortran_file_open(file_id, file_name, file_status, err,error,*999)
        section_start=.false.
        file_end=.false.

        DO WHILE(.NOT.file_end)
          CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)

          !check the beginning of field section in exnode files
          IF((.NOT.section_start).AND.(verify(cmiss_keyword_fields,line)==0)) THEN
            section_start=.true.
          ENDIF

          !check whether it is a new header for another group of elements
          IF(section_start.AND.(verify(cmiss_keyword_fields,line)==0)) THEN

            !collect header information
            pos=index(line,cmiss_keyword_fields)
            line=remove(line,1, pos+len_trim(cmiss_keyword_fields)-1)
            idx_field=string_to_integer(line, err,error)
            IF(idx_field/=number_of_fields) CALL flagerror("find different field number in exnode files",err,error,*999)
            idx_comp=0
            DO idx_field=1,number_of_fields
              CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
              IF(idx_exnode==0) THEN
                list_str(idx_field)=line
                pos=index(line,cmiss_keyword_components)
                line=remove(line, 1, pos+len_trim(cmiss_keyword_components)-1)
                number_of_comps=string_to_integer(line, err,error)
                total_number_of_comps=total_number_of_comps+number_of_comps
              ELSE
             IF(list_str(idx_field)/=line) CALL flagerror("find different field information in exnode files", &
              & err,error,*999)
              ENDIF
              pos=index(line,cmiss_keyword_components)
              line=remove(line,1, pos+len_trim(cmiss_keyword_components)-1)
              number_of_comps=string_to_integer(line, err,error)
              DO idx_comp=1,number_of_comps
                CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
              ENDDO !idx_comp1
            ENDDO !idx_field
          ENDIF !START_OF_FIELD_SECTION==.TRUE..AND.(VERIFY(CMISS_KEYWORD_FIELD,LINE)==0)
        ENDDO !(FILE_END==.FALSE.)
        CALL field_io_fortran_file_close(file_id, err,error,*999)
        idx_exnode=idx_exnode+1
      ENDDO !idx_exnode<NUMBER_OF_EXNODE_FILES
    ENDIF !MASTER_COMPUTATIONAL_NUMBER==my_computational_node_number

    idx_comp1=0
    DO idx_field=1,number_of_fields
      IF(ASSOCIATED(field)) NULLIFY(field)
      !Start to create a default (geometric) field on the region
      CALL field_create_start(idx_field,region,field,err,error,*999)
      !always has one field variable in one field during reading
      CALL field_number_of_variables_set(field,1,err,error,*999)
      !Set the decomposition to use
      CALL field_mesh_decomposition_set(field,decomposition,err,error,*999)
      !Set the number of components for this field
      CALL field_number_of_components_set(field,field_u_variable_type,components_in_fields(idx_field),err,error,*999)
      DO idx_comp=1, components_in_fields(idx_field)
        idx_comp1=idx_comp1+1
        !Set the domain to be used by the field components
        CALL field_component_mesh_component_set(field,1,idx_comp,mesh_components_of_field_components(idx_comp1),err,error,*999)
      ENDDO
      !Set the scaling factor
      CALL field_scaling_type_set(field, field_scaling_type, err, error, *999)

      IF(master_computational_number==my_computational_node_number) THEN
        CALL field_io_field_info(list_str(idx_field), field_io_field_label, fieldtype, err, error, *999)
      ENDIF
      CALL mpi_bcast(fieldtype,1,mpi_logical,master_computational_number,mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
      !Set FIELD TYPE
      CALL field_type_set(field, fieldtype, err, error, *999)
      !Finish creating the field
      CALL field_create_finish(field,err,error,*999)
    ENDDO

    IF(master_computational_number==my_computational_node_number) THEN
      CALL checked_deallocate( list_str )
    ENDIF

    file_end=.true.
    idx_exnode=-1
    file_id=1030
    file_status="OLD"

    !broadcasting total_number_of_comps
    CALL mpi_bcast(total_number_of_comps,1,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
    CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)

    CALL reallocate( list_dev_pos, total_number_of_comps, &
      & "Could not allocate memory for nodal derivative position in field components", err, error, *999 )

    DO WHILE(idx_exnode<number_of_exnode_files)

      CALL mpi_bcast(file_end,1,mpi_logical,master_computational_number,mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)

      IF(file_end) THEN
        idx_exnode=idx_exnode+1
        INQUIRE(unit=file_id, opened=file_open)
        IF(file_open) CALL field_io_fortran_file_close(file_id, err,error,*999)
        IF(idx_exnode>=number_of_exnode_files) EXIT
      ENDIF

      !IF(MASTER_COMPUTATIONAL_NUMBER/=my_computational_node_number) PRINT * , idx_exnode

      !goto the start of mesh part
      IF(master_computational_number==my_computational_node_number) THEN

        IF(file_end) THEN
          file_id=1030+idx_exnode
          !checking the next file
          file_name=name//".part"//trim(number_to_vstring(idx_exnode,"*",err,error))//".exnode"
          CALL field_io_fortran_file_open(file_id, file_name, file_status, err,error,*999)
          file_end=.false.
          section_start=.false.
          node_section=.false.
          !idx_exnode=idx_exnode+1
        ENDIF

        IF((.NOT.file_end).AND.(.NOT.section_start))  THEN
          !find a new header
          DO WHILE(verify(cmiss_keyword_fields,line)/=0)
            CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
          ENDDO
          section_start=.true.
        ENDIF

        !have not touched the end
        IF((.NOT.file_end).AND.section_start.AND.(.NOT.node_section)) THEN
          pos=index(line,cmiss_keyword_fields)
          line=remove(line,1, pos+len_trim(cmiss_keyword_fields)-1)
          !number_of_fields=STRING_TO_INTEGER(LINE, ERR, ERROR)
          total_number_of_devs=0
          idx_comp1=0
          idx_dev1=0
          DO idx_field=1, number_of_fields
            CALL field_io_fortran_file_read_string(file_id, line, file_end, err, error,*999)
            pos=index(line,cmiss_keyword_components)
            line=remove(line,1, pos+len_trim(cmiss_keyword_components)-1)
            number_of_comps=string_to_integer(line, err, error)
            !total_number_of_comps=total_number_of_comps+number_of_comps

            DO idx_comp=1, number_of_comps
              CALL field_io_fortran_file_read_string(file_id, line, file_end, err, error,*999)
              pos=index(line,".")
              line=remove(line,1, pos+1)
              pos=index(line,cmiss_keyword_value_index)
              line=remove(line,1, pos+len_trim(cmiss_keyword_value_index)-1)
              pos=index(line,",")
              line1=extract(line,1,pos-1)
              idx_comp1=idx_comp1+1
              list_dev_pos(idx_comp1)=string_to_integer(line1, err, error)

              pos=index(line,cmiss_keyword_derivative)
              line=remove(line,1, pos+len_trim(cmiss_keyword_derivative)-1)
              pos=index(line,"(")
              line1=extract(line,1,pos-1)
              number_of_devs=string_to_integer(line, err, error)+1
              total_number_of_devs=total_number_of_devs+number_of_devs

              IF(ALLOCATED(list_dev)) THEN
                CALL reallocate( tmp_pointer, total_number_of_devs-number_of_devs, &
                  & "Could not allocate temporary memory for nodal derivative index in master node", err,error,*999)
                tmp_pointer(:)=list_dev(:)

                CALL reallocate( list_dev, total_number_of_devs, &
                  & "Could not allocate temporary memory for nodal derivative index in master node", err,error,*999)
                list_dev(1:total_number_of_devs-number_of_devs)=tmp_pointer(:)

                DEALLOCATE(tmp_pointer)
              ELSE
                CALL reallocate( list_dev, total_number_of_devs, &
                  & "Could not allocate memory for nodal derivative index", err, error, *999)
              ENDIF

              !print *, idx_dev1, NO_PART_DERIV, LIST_DEV

              IF(number_of_devs<=1) THEN
                idx_dev1=idx_dev1+1
                list_dev(idx_dev1)=no_part_deriv
                !print *, idx_dev1, NO_PART_DERIV, LIST_DEV
              ELSE
                pos=index(line,"(")
                line=remove(line,1, pos)
                pos=index(line,")")
                line=remove(line,pos, len(line))
                idx_dev1=idx_dev1+1
                list_dev(idx_dev1)=no_part_deriv
                DO idx_dev=2, number_of_devs-1
                  idx_dev1=idx_dev1+1
                  pos=index(line,",")
                  line1=extract(line, 1, pos-1)
                  line=remove(line, 1, pos)
                  list_dev(idx_dev1)=field_io_derivative_info(line1, err,error)
                ENDDO
                idx_dev1=idx_dev1+1
                list_dev(idx_dev1)=field_io_derivative_info(line, err,error)
              ENDIF
            ENDDO !idx_comp
            CALL field_io_fortran_file_read_string(file_id, line, file_end, err, error,*999)
            node_section=.true.
          ENDDO !idx_field
        ENDIF  !FILE_END==.FALSE..AND.SECTION_START=.TRUE..AND.NODE_SECTION=.FALSE.
      ENDIF !MASTER_COMPUTATIONAL_NUMBER

      !broadcasting total_number_of_devs
      CALL mpi_bcast(total_number_of_devs,1,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)

      IF(master_computational_number/=my_computational_node_number) THEN
        CALL reallocate( list_dev, total_number_of_devs, &
          & "Could not allocate memory for nodal derivative index in non-master node", err, error, *999 )
      ENDIF

      CALL reallocate( list_dev_value, total_number_of_devs, &
        & "Could not allocate memory for nodal derivative index in non-master node", err, error, *999 )

      !broadcasting total_number_of_comps
      CALL mpi_bcast(list_dev_pos,total_number_of_comps,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
      !broadcasting total_number_of_devs
      CALL mpi_bcast(list_dev,total_number_of_devs,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)

      !goto the start of mesh part
      IF(master_computational_number==my_computational_node_number) THEN

        !have not touched the end
        IF((.NOT.file_end).AND.section_start.AND.node_section) THEN

          IF(verify(cmiss_keyword_node, line)==0) THEN
            pos=index(line,cmiss_keyword_node)
            line=remove(line,1, pos+len_trim(cmiss_keyword_node)-1)
            nodal_user_number=string_to_integer(line, err, error)
            idx_comp1=1
            DO idx_comp=1, number_of_comps-1
              IF(list_dev_pos(idx_comp+1)-list_dev_pos(idx_comp)<=number_nodal_value_lines) THEN
                CALL field_io_fortran_file_read_string(file_id, line, file_end, err, error,*999)
                CALL string_to_muti_reals_vs(line, list_dev_pos(idx_comp+1)-list_dev_pos(idx_comp), list_dev_value, &
                  & list_dev_pos(idx_comp), err,error, *999)
              ELSE
                number_of_lines=(list_dev_pos(idx_comp+1)-list_dev_pos(idx_comp))/number_nodal_value_lines
                DO idx_nodal_line=1, number_of_lines
                  CALL field_io_fortran_file_read_string(file_id, line, file_end, err, error,*999)
                  CALL string_to_muti_reals_vs(line, number_nodal_value_lines, list_dev_value, list_dev_pos(idx_comp)+ &
                    & (idx_nodal_line-1)*number_nodal_value_lines, err,error, *999)
                ENDDO
                CALL field_io_fortran_file_read_string(file_id, line, file_end, err, error,*999)
                CALL string_to_muti_reals_vs(line, (list_dev_pos(idx_comp+1)-list_dev_pos(idx_comp))- &
                  & number_nodal_value_lines*number_of_lines, list_dev_value, list_dev_pos(idx_comp)+ &
                  & (idx_nodal_line-1)*number_nodal_value_lines, err,error, *999)
              ENDIF
            ENDDO
            !IF((total_number_of_devs-LIST_DEV_POS(idx_comp)+1)<=NUMBER_NODAL_VALUE_LINES) THEN
            !  CALL FIELD_IO_FORTRAN_FILE_READ_STRING(FILE_ID, LINE, FILE_END, ERR, ERROR,*999)
            !  CALL STRING_TO_MUTI_REALS_VS(LINE, total_number_of_devs-LIST_DEV_POS(idx_comp)+1, LIST_DEV_VALUE, LIST_DEV_POS(idx_comp), ERR,ERROR, *999)
            !ELSE
            !  NUMBER_OF_LINES=(total_number_of_devs-LIST_DEV_POS(idx_comp)+1)/NUMBER_NODAL_VALUE_LINES
            !  DO idx_nodal_line=1, NUMBER_OF_LINES
            !    CALL FIELD_IO_FORTRAN_FILE_READ_STRING(FILE_ID, LINE, FILE_END, ERR, ERROR,*999)
            !    CALL STRING_TO_MUTI_REALS_VS(LINE, NUMBER_NODAL_VALUE_LINES, LIST_DEV_VALUE, LIST_DEV_POS(idx_comp)+(idx_nodal_line-1)*NUMBER_NODAL_VALUE_LINES, ERR,ERROR, *999)
            !  ENDDO
            !  CALL FIELD_IO_FORTRAN_FILE_READ_STRING(FILE_ID, LINE, FILE_END, ERR, ERROR,*999)
            !  CALL STRING_TO_MUTI_REALS_VS(LINE, (LIST_DEV_POS(idx_comp+1)-LIST_DEV_POS(idx_comp))-NUMBER_NODAL_VALUE_LINES*NUMBER_OF_LINES, LIST_DEV_VALUE, LIST_DEV_POS(idx_comp)+(idx_nodal_line-1)*NUMBER_NODAL_VALUE_LINES, ERR,ERROR, *999)
            !ENDIF
            CALL field_io_fortran_file_read_string(file_id, line, file_end, err, error,*999)
            CALL string_to_muti_reals_vs(line, total_number_of_devs-list_dev_pos(idx_comp)+1, list_dev_value, &
              & list_dev_pos(idx_comp), err,error, *999)
          ELSE
            CALL flagerror("The position of nodal information in exenode files is not correct",err, error,*999)
            node_section=.false.
          ENDIF !(VERIFY(CMISS_KEYWORD_NODE , LINE)==0)
          IF(.NOT.file_end) THEN
            CALL field_io_fortran_file_read_string(file_id, line, file_end, err, error,*999)
            IF(verify(cmiss_keyword_node, line)/=0) node_section=.false.
          ENDIF
        ENDIF !FILE_END==.FALSE..AND.SECTION_START=.TRUE..AND.NODE_SECTION=.TRUE.
      ENDIF  !(MASTER_COMPUTATIONAL_NUMBER==my_computational_node_number)

      !broadcasting total_number_of_devs
      CALL mpi_bcast(list_dev_value,total_number_of_devs,mpi_real8,master_computational_number,mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
      CALL mpi_bcast(nodal_user_number,1,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)

      !IF(MASTER_COMPUTATIONAL_NUMBER/=my_computational_node_number) THEN
      print *, "user number:"
      print *, nodal_user_number
      print *, list_dev_value
      !ENDIF


      idx_comp1=0
      idx_dev1=0
      idx_variable=1
      DO idx_field=1,number_of_fields
        IF(ASSOCIATED(field)) NULLIFY(field)
        field=>region%FIELDS%FIELDS(idx_field)%PTR
        DO idx_comp=1, components_in_fields(idx_field)
          idx_comp1=idx_comp1+1
          domain_nodes=>field%VARIABLES(idx_variable)%COMPONENTS(idx_comp)%DOMAIN%TOPOLOGY%NODES
          node_in_domain=.false.
          DO idx_node=1,domain_nodes%NUMBER_OF_NODES
            !IF(DOMAIN_NODES%NODES(idx_node)%GLOBAL_NUMBER==USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER(NODAL_USER_NUMBER)) THEN
            IF(domain_nodes%NODES(idx_node)%USER_NUMBER==nodal_user_number) THEN
              node_in_domain=.true.
              nodal_local_number=idx_node
            ENDIF
          ENDDO

          IF(node_in_domain) THEN
            IF(idx_comp1>=number_of_components) THEN
              DO idx_dev=1, total_number_of_devs-list_dev_pos(idx_comp1)+1
                idx_dev1=idx_dev1+1
                !Set the domain to be used by the field components
                !Default to version 1 of each node derivative
                CALL field_parameter_set_update_node(field,field_values_set_type,1, list_dev(idx_dev1), &
                    &nodal_local_number, idx_comp, idx_variable, list_dev_value(idx_dev1),&
                    &err, error, *999)
                !print *, "n--n"
              ENDDO !idx_dev
            ELSE
              DO idx_dev=1, list_dev_pos(idx_comp1+1)-list_dev_pos(idx_comp1)
                idx_dev1=idx_dev1+1
                !Set the domain to be used by the field components
                !Default to version 1 of each node derivative
                CALL field_parameter_set_update_node(field,field_values_set_type,1, list_dev(idx_dev1), &
                    &nodal_local_number, idx_comp, idx_variable, list_dev_value(idx_dev1),&
                    &err, error, *999)
                !print *, "n--n"
              ENDDO !idx_dev
            ENDIF  !idx_comp1
          ENDIF !NODE_IN_DOMAIN
        ENDDO !idx_comp
      ENDDO !idx_field
    ENDDO !idx_exnode<NUMBER_OF_EXELEM_FILES

    !print *, "out of loop"

    DO idx_field=1,number_of_fields
       IF(ASSOCIATED(field)) NULLIFY(field)
       field=>region%FIELDS%FIELDS(idx_field)%PTR
       DO variable_idx=1,field%NUMBER_OF_VARIABLES
         variable_type=field%VARIABLES(variable_idx)%VARIABLE_TYPE
         CALL field_parameter_set_update_start(field,variable_type,field_values_set_type,err,error,*999)
         CALL field_parameter_set_update_finish(field,variable_type,field_values_set_type,err,error,*999)
       ENDDO !variable_idx
    ENDDO

    IF(ALLOCATED(tmp_pointer)) DEALLOCATE(tmp_pointer)
    IF(ALLOCATED(list_dev_value)) DEALLOCATE(list_dev_value)
    IF(ALLOCATED(list_dev)) DEALLOCATE(list_dev)
    IF(ALLOCATED(list_dev_pos)) DEALLOCATE(list_dev_pos)
    IF(ALLOCATED(list_str)) DEALLOCATE(list_str)

    exits("FIELD_IO_CREATE_FIELDS")
    RETURN
999 errorsexits("FIELD_IO_CREATE_FIELDS",err,error)
    RETURN 1
  END SUBROUTINE field_io_create_fields


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

  SUBROUTINE field_io_create_decompistion(DECOMPOSITION, DECOMPOSITION_USER_NUMBER, DECOMPOSITION_METHOD, MESH, &
    & number_of_domains, err, error, *)
    !Argument variables
    TYPE(decomposition_type), POINTER :: DECOMPOSITION
    INTEGER(INTG), INTENT(IN) :: DECOMPOSITION_USER_NUMBER
    INTEGER(INTG), INTENT(IN) :: DECOMPOSITION_METHOD
    TYPE(mesh_type), POINTER :: MESH
    INTEGER(INTG), INTENT(IN) :: NUMBER_OF_DOMAINS
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

    IF(.NOT.ASSOCIATED(mesh)) THEN
      CALL flagerror("mesh is NOT associated before decomposing the mesh",err,error,*999)
      GOTO 999
    ENDIF

    enters("FIELD_IO_CREATE_DECOMPISTION",err,error,*999)
    !Create a decomposition
    CALL decomposition_create_start(decomposition_user_number,mesh,decomposition,err,error,*999)
    !Set the decomposition to be a general decomposition with the specified number of domains
    CALL decomposition_type_set(decomposition,decomposition_method,err,error,*999)
    CALL decomposition_number_of_domains_set(decomposition,number_of_domains,err,error,*999)
    CALL decomposition_create_finish(decomposition,err,error,*999)

    exits("FIELD_IO_CREATE_DECOMPISTION")
    RETURN
999 errorsexits("FIELD_IO_CREATE_DECOMPISTION",err,error)
    RETURN 1
  END SUBROUTINE field_io_create_decompistion


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

  SUBROUTINE field_io_fields_import(NAME, METHOD, REGION, MESH, MESH_USER_NUMBER, DECOMPOSITION, DECOMPOSITION_USER_NUMBER, &
    &decomposition_method, field_values_set_type, field_scaling_type, err, error, *)
    !Argument variables
    TYPE(varying_string), INTENT(IN) :: NAME
    TYPE(varying_string), INTENT(IN) :: METHOD
    TYPE(region_type), POINTER :: REGION
    TYPE(mesh_type), POINTER :: MESH
    INTEGER(INTG), INTENT(IN) :: MESH_USER_NUMBER
    TYPE(decomposition_type), POINTER :: DECOMPOSITION
    INTEGER(INTG), INTENT(IN) :: DECOMPOSITION_USER_NUMBER
    INTEGER(INTG), INTENT(IN) :: DECOMPOSITION_METHOD
    INTEGER(INTG), INTENT(IN) :: FIELD_VALUES_SET_TYPE
    INTEGER(INTG), INTENT(IN) :: FIELD_SCALING_TYPE
    !TYPE(BASIS_FUNCTIONS_TYPE), POINTER :: BASES !< bases function
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: my_computational_node_number !local process number
    INTEGER(INTG) :: computational_node_numbers   !total process numbers
    INTEGER(INTG) :: MASTER_COMPUTATIONAL_NUMBER  !master computational number
    INTEGER(INTG) :: NUMBER_OF_FIELDS
    INTEGER(INTG) :: NUMBER_OF_EXNODE_FILES
    !INTEGER(INTG), ALLOCATABLE :: USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER(:)
    INTEGER(INTG), ALLOCATABLE :: MESH_COMPONENTS_OF_FIELD_COMPONENTS(:)
    INTEGER(INTG), ALLOCATABLE :: COMPONENTS_IN_FIELDS(:)

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

    !Get the number of computational nodes
    computational_node_numbers=computational_nodes_number_get(err,error)
    IF(err/=0) GOTO 999
    !Get my computational node number
    my_computational_node_number=computational_node_number_get(err,error)
    IF(err/=0) GOTO 999

    master_computational_number=0

    IF(method=="FORTRAN") THEN
      CALL field_io_import_global_mesh(name, region, mesh, mesh_user_number, master_computational_number, &
          & my_computational_node_number, &!USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER,
          &mesh_components_of_field_components, &
          & components_in_fields, number_of_fields, number_of_exnode_files, err, error, *999)

      CALL field_io_create_decompistion(decomposition, decomposition_user_number, decomposition_method, mesh, &
          &computational_node_numbers, err, error, *999)

      CALL field_io_create_fields(name, region, decomposition, field_values_set_type, number_of_fields, &
          !&USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER,
          &mesh_components_of_field_components, components_in_fields, &
          & number_of_exnode_files, master_computational_number, my_computational_node_number, field_scaling_type, &
          & err, error, *999)
    ELSE IF(method=="MPIIO") THEN
      CALL flagerror("MPI IO has not been implemented",err,error,*999)
    ELSE
      CALL flagerror("Unknown method!",err,error,*999)
    ENDIF

    !IF(ALLOCATED(USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER)) DEALLOCATE(USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER)
    CALL checked_deallocate( mesh_components_of_field_components )
    CALL checked_deallocate( components_in_fields )
    !IF(ALLOCATED(LIST_FIELD_TYPE)) DEALLOCATE(LIST_FIELD_TYPE)

    exits("FIELD_IO_FIELDS_IMPORT")
    RETURN
999 errorsexits("FIELD_IO_FIELDS_IMPORT",err,error)
    RETURN 1
  END SUBROUTINE field_io_fields_import

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

  SUBROUTINE field_io_fill_basis_info(INTERPOLATION_XI, LIST_STR, NUMBER_OF_COMPONENTS, ERR, ERROR, *)
    !Argument variables
    INTEGER(INTG), INTENT(INOUT) :: INTERPOLATION_XI(:,:)
    TYPE(varying_string), INTENT(INOUT) :: LIST_STR(:)
    INTEGER(INTG), INTENT(IN) :: NUMBER_OF_COMPONENTS ! number of components
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LINE, LINE1
    INTEGER(INTG) :: idx_comp, pos
    INTEGER(INTG) :: num_interp, INTERPOLATION_TYPE

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

    DO idx_comp=1,number_of_components
      num_interp=0
      line=list_str(idx_comp)
      DO WHILE(verify("*",line)==0)
        num_interp=num_interp+1
        pos=index(line,"*")
        line1=extract(line, 1, pos)
        line=remove(line,1,pos)
        CALL fieldio_translatelabelintointerpolationtype(interpolation_type, line, err, error, *999)
        interpolation_xi(idx_comp, num_interp)=interpolation_type
      ENDDO
      num_interp=num_interp+1
      line1=extract(line, 1, pos)
      line=remove(line,1,pos)
      CALL fieldio_translatelabelintointerpolationtype(interpolation_type, line, err, error, *999)
      interpolation_xi(idx_comp, num_interp)=interpolation_type
    ENDDO

    exits("FIELD_IO_FILL_BASIS_INFO")
    RETURN
999 errorsexits("FIELD_IO_FILL_BASIS_INFO",err,error)
    RETURN 1
  END SUBROUTINE field_io_fill_basis_info


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

  SUBROUTINE field_io_import_global_mesh(NAME, REGION, MESH, MESH_USER_NUMBER, MASTER_COMPUTATIONAL_NUMBER, &
    & my_computational_node_number, &!USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER,
    &mesh_components_of_field_components, &
    & components_in_fields, number_of_fields, number_of_exnode_files, err, error, *)
    !Argument variables
    TYPE(varying_string), INTENT(IN):: NAME
    TYPE(mesh_type), POINTER :: MESH
    TYPE(region_type), POINTER :: REGION
    INTEGER(INTG), INTENT(IN) :: MESH_USER_NUMBER
    INTEGER(INTG), INTENT(IN) :: MASTER_COMPUTATIONAL_NUMBER
    INTEGER(INTG), INTENT(IN) :: my_computational_node_number
    !INTEGER(INTG), INTENT(INOUT), ALLOCATABLE :: USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER(:)
    INTEGER(INTG), INTENT(INOUT), ALLOCATABLE :: MESH_COMPONENTS_OF_FIELD_COMPONENTS(:)
    INTEGER(INTG), INTENT(INOUT), ALLOCATABLE :: COMPONENTS_IN_FIELDS(:)
    INTEGER(INTG), INTENT(INOUT) :: NUMBER_OF_FIELDS
    INTEGER(INTG), INTENT(INOUT) :: NUMBER_OF_EXNODE_FILES
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string), ALLOCATABLE :: LIST_STR(:)
    TYPE(basis_type), POINTER :: BASIS
    TYPE(nodes_type), POINTER :: NODES
    TYPE(mesh_elements_type_ptr_type), ALLOCATABLE :: ELEMENTS_PTR(:)
    TYPE(varying_string) :: FILE_NAME, FILE_STATUS, LINE
    TYPE(varying_string) :: CMISS_KEYWORD_FIELDS, CMISS_KEYWORD_ELEMENT, CMISS_KEYWORD_NODE, CMISS_KEYWORD_COMPONENTS
    TYPE(varying_string) :: CMISS_KEYWORD_SHAPE, CMISS_KEYWORD_SCALE_FACTOR_SETS, CMISS_KEYWORD_NODES, CMISS_KEYWORD_SCALE_FACTORS
    INTEGER(INTG), PARAMETER :: NUMBER_NODAL_LINES=3, number_scaling_factors_in_line=5
    INTEGER(INTG), ALLOCATABLE :: LIST_ELEMENT_NUMBER(:), LIST_ELEMENTAL_NODES(:), LIST_COMP_NODAL_INDEX(:,:)
    INTEGER(INTG), ALLOCATABLE :: MESH_COMPONENT_LOOKUP(:,:), INTERPOLATION_XI(:,:), LIST_COMP_NODES(:)!LIST_FIELD_COMPONENTS(:),
    INTEGER(INTG), ALLOCATABLE :: USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER(:)
    INTEGER(INTG) :: FILE_ID, NUMBER_OF_EXELEM_FILES, NUMBER_OF_ELEMENTS, NUMBER_OF_NODES, NUMBER_OF_DIMENSIONS
    INTEGER(INTG) :: NUMBER_OF_MESH_COMPONENTS, NUMBER_OF_COMPONENTS, NUMBER_SCALING_FACTOR_LINES
    INTEGER(INTG) :: GLOBAL_ELEMENT_NUMBER
    INTEGER(INTG) :: MPI_IERROR
    INTEGER(INTG) :: SHAPE_INDEX(shape_size)
    INTEGER(INTG) :: idx_comp, idx_comp1, pos, idx_node, idx_node1, idx_field, idx_elem, idx_exnode, idx_exelem, number_of_comp
    INTEGER(INTG) :: idx_basis, number_of_node, number_of_scalesets, idx_scl, idx_mesh_comp, current_mesh_comp, num_scl,&
      & num_scl_line
    LOGICAL :: FILE_EXIST, START_OF_ELEMENT_SECTION, FIELD_SECTION, SECTION_START, FILE_END, FILE_OPEN

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

    !checking region pointer
    IF(.NOT.ASSOCIATED(region)) THEN
      CALL flagerror("region is not associated",err,error,*999)
      GOTO 999
    ENDIF

    !checking mesh pointer
    IF(ASSOCIATED(mesh)) THEN
      CALL flagerror("mesh is associated, pls release the memory first",err,error,*999)
      GOTO 999
    ENDIF

    IF(.NOT.region%REGION_FINISHED) THEN
      CALL flagerror("region is not finished",err,error,*999)
      GOTO 999
    ENDIF

    !IF(ASSOCIATED(BASES%BASES)) THEN
    !   CALL FlagError("bases are associated, pls release the memory first",ERR,ERROR,*999)
    !   GOTO 999
    !ENDIF
    !BASES%NUMBER_BASIS_FUNCTIONS=0
    number_of_dimensions=region%COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS
    file_status="OLD"
    cmiss_keyword_shape="Shape.  Dimension="//trim(number_to_vstring(number_of_dimensions,"*",err,error))
    cmiss_keyword_element="Element:"
    cmiss_keyword_components="#Components="
    cmiss_keyword_node="Node:"
    cmiss_keyword_nodes="#Nodes="
    cmiss_keyword_fields="#Fields="
    cmiss_keyword_scale_factor_sets="#Scale factor sets="
    cmiss_keyword_scale_factors="#Scale factors="

    CALL mesh_create_start(mesh_user_number,region,number_of_dimensions,mesh,err,error,*999)

    !calculate the number of elements, number of fields and number of field components
    IF(master_computational_number==my_computational_node_number) THEN

      !the file name has to start from zero in an ascended order without break
      idx_exelem=0
      idx_elem=0
      number_of_components=0
      file_name=name//".part"//trim(number_to_vstring(idx_exelem,"*",err,error))//".exelem"
      INQUIRE(file=char(file_name), exist=file_exist)
      IF(.NOT.file_exist) THEN
        CALL flagerror("exelem files can be found, pls check again",err,error,*999)
        !GOTO 999
      ENDIF
      DO WHILE(file_exist)

        file_id=1030+idx_exelem
        CALL field_io_fortran_file_open(file_id, file_name, file_status, err,error,*999)
        start_of_element_section=.false.
        field_section=.false.

        CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error,*999)
        DO WHILE(.NOT.file_end)

          !check the beginning of element section
          IF((.NOT.start_of_element_section).AND.(verify(cmiss_keyword_shape,line)==0)) THEN
            start_of_element_section=.true.
          ENDIF

          !count how many elements
          IF(start_of_element_section.AND.(verify(cmiss_keyword_element,line)==0)) idx_elem=idx_elem+1

          !check whether they have same numbers of fields
          IF(start_of_element_section.AND.(verify(cmiss_keyword_fields,line)==0)) THEN
            idx_field=0
            idx_comp=0
            field_section=.true.
            pos=index(line,cmiss_keyword_fields)
            line=remove(line,1, pos+len_trim(cmiss_keyword_fields)-1)
            IF(idx_exelem==0) THEN
              number_of_fields=string_to_integer(line, err,error)
            ELSE
              IF(number_of_fields/=string_to_integer(line, err,error)) THEN
                CALL flagerror("find different number of fields in the exelem files",err,error,*999)
                !GOTO 999
              ENDIF
            ENDIF

            IF(.NOT.ALLOCATED(components_in_fields)) THEN
              CALL reallocate( components_in_fields, number_of_fields, &
                & "can not allocate the memory for outputing components in field", err, error, *999 )
            ENDIF
          ENDIF !START_OF_ELEMENT_SECTION==.TRUE..AND.VERIFY(CMISS_KEYWORD_FIELDS,LINE)==0

          !check whether they have same numbers of field components
          IF(field_section.AND.start_of_element_section.AND.(verify(cmiss_keyword_components,line)==0)) THEN
            idx_field=idx_field+1
            pos=index(line,cmiss_keyword_components)
            line=remove(line,1, pos+len_trim(cmiss_keyword_components)-1)
            idx_comp1=string_to_integer(line, err,error)
            idx_comp=idx_comp+idx_comp1
            IF(idx_field>=number_of_fields) THEN
              IF(idx_exelem==0) THEN
                number_of_components=idx_comp
                components_in_fields(idx_field)=idx_comp
              ELSE
                IF(number_of_components/=idx_comp) THEN
                  CALL flagerror("find different total number of components in the exelem files",err,error,*999)
                  !GOTO 999
                ENDIF
              ENDIF
              field_section=.false.
            ENDIF
            IF(idx_exelem==0) THEN
              components_in_fields(idx_field)=idx_comp1
            ELSE
              IF(components_in_fields(idx_field)/=idx_comp1) THEN
                CALL flagerror("find different number of components in one field in the exelem files",err,error,*999)
                !GOTO 999
              ENDIF
            ENDIF
          ENDIF !FIELD_SECTION==.TRUE..AND.START_OF_ELEMENT_SECTION==.TRUE..AND.VERIFY(CMISS_KEYWORD_COMPONENTS,LINE
          CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error,*999)
        ENDDO !(FILE_END==.FALSE.)

        CALL field_io_fortran_file_close(file_id, err,error,*999)
        !checking the next file
        idx_exelem=idx_exelem+1
        file_name=name//".part"//trim(number_to_vstring(idx_exelem,"*",err,error))//".exelem"
        INQUIRE(file=char(file_name), exist=file_exist)
      ENDDO !FILE_EXIST==.TRUE.
      !CALL WRITE_STRING_VALUE(GENERAL_OUTPUT_TYPE,"  Total number of exelment files = ",idx_exelem, ERR,ERROR,*999)
      number_of_elements=idx_elem
      number_of_exelem_files=idx_exelem
    ENDIF !MASTER_COMPUTATIONAL_NUMBER==my_computational_node_number

    !broadcasting the number of components in each field
    CALL mpi_bcast(number_of_fields,1,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
    CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
    IF(master_computational_number/=my_computational_node_number) THEN
      CALL reallocate( components_in_fields, number_of_fields, &
          & "can not allocate the memory for outputing components in field", err, error, *999 )
      !IF(ALLOCATED(LIST_FIELD_TYPE)) DEALLOCATE(LIST_FIELD_TYPE)
      !ALLOCATE(LIST_FIELD_TYPE(NUMBER_OF_FIELDS), STAT=ERR)
      !IF(ERR/=0) CALL FlagError("can not allocate the memory for list of field types",ERR,ERROR,*999)
      !LIST_FIELD_TYPE(:)=0
    ENDIF
    CALL mpi_bcast(components_in_fields,number_of_fields,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
    CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
    !CALL MPI_BCAST(LIST_FIELD_TYPE,NUMBER_OF_FIELDS,MPI_INTEGER,MASTER_COMPUTATIONAL_NUMBER,MPI_COMM_WORLD,MPI_IERROR)
    !CALL MPI_ERROR_CHECK("MPI_BCAST",MPI_IERROR,ERR,ERROR,*999)
    !broadcasting the number of elements
    CALL mpi_bcast(number_of_elements,1,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
    CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
    CALL mesh_number_of_elements_set(mesh,number_of_elements,err,error,*999)

    !calculate the number of nodes
    IF(master_computational_number==my_computational_node_number) THEN
      !the file name has to start from zero in a ascended order without break
      idx_exnode=0
      idx_node=0
      file_name=name//".part"//trim(number_to_vstring(idx_exnode,"*",err,error))//".exnode"
      INQUIRE(file=char(file_name), exist=file_exist)
      IF(.NOT.file_exist) THEN
        CALL flagerror("exnode files can be found, pls check again",err,error,*999)
        !GOTO 999
      ENDIF
      DO WHILE(file_exist)
        file_id=1030+idx_exnode
        CALL field_io_fortran_file_open(file_id, file_name, file_status, err,error,*999)
        file_end=.false.

        DO WHILE(.NOT.file_end)
          CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
          IF((.NOT.file_end).AND.verify(cmiss_keyword_node,line)==0) idx_node=idx_node+1
        ENDDO !(FILE_END==.FALSE.)

        CALL field_io_fortran_file_close(file_id, err,error,*999)
        !checking the next file
        idx_exnode=idx_exnode+1
        file_name=name//".part"//trim(number_to_vstring(idx_exnode,"*",err,error))//".exnode"
        INQUIRE(file=char(file_name), exist=file_exist)
      ENDDO !FILE_EXIST==.TRUE.
      !CALL WRITE_STRING_VALUE(GENERAL_OUTPUT_TYPE,"  Total number of exnode files = ",idx_exnode, ERR,ERROR,*999)
      number_of_nodes=idx_node
      number_of_exnode_files=idx_exnode
    ENDIF !MASTER_COMPUTATIONAL_NUMBER==my_computational_node_number

    CALL mpi_bcast(number_of_exnode_files,1,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
    CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
    !broadcasting the number of nodes
    CALL mpi_bcast(number_of_nodes,1,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
    CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
    NULLIFY(nodes)
    CALL nodes_create_start(region,number_of_nodes,nodes,err,error,*999)

    !collect the nodal numberings (nodal labels) to change the nodal user number by reading exnode files
    CALL reallocate( user_nodal_number_map_global_nodal_number, number_of_nodes, &
      & "can not allocate list of nodal number.", err, error, *999 )
    IF(master_computational_number==my_computational_node_number) THEN
      !the file name has to start from zero in a ascended order without break
      idx_node=1
      DO idx_exnode=0, number_of_exnode_files-1
        file_id=1030+idx_exnode
        file_name=name//".part"//trim(number_to_vstring(idx_exnode,"*",err,error))//".exnode"
        CALL field_io_fortran_file_open(file_id, file_name, file_status, err,error,*999)
        file_end=.false.
        DO WHILE(.NOT.file_end)
          CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
          IF((.NOT.file_end).AND.verify(cmiss_keyword_node,line)==0) THEN
            pos=index(line,cmiss_keyword_node)
            line=remove(line,1, pos+len_trim(cmiss_keyword_node)-1)
            user_nodal_number_map_global_nodal_number(idx_node)=string_to_integer(line, err, error)
            idx_node=idx_node+1
          ENDIF !VERIFY(CMISS_KEYWORD,LINE)==0
        ENDDO !(FILE_END==.FALSE.)
        CALL field_io_fortran_file_close(file_id, err,error,*999)
      ENDDO !FILE_EXIST==.TRUE.
      CALL list_sort(user_nodal_number_map_global_nodal_number, err, error, *999)
    ENDIF !MASTER_COMPUTATIONAL_NUMBER==my_computational_node_number

    !broadcast the nodal numberings (nodal labels)
    CALL mpi_bcast(user_nodal_number_map_global_nodal_number,number_of_nodes,mpi_integer,master_computational_number, &
        & mpi_comm_world,mpi_ierror)
    CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
    DO idx_node=1, number_of_nodes
      IF(idx_node/=user_nodal_number_map_global_nodal_number(idx_node)) CALL nodes_user_number_set(nodes,idx_node, &
          & user_nodal_number_map_global_nodal_number(idx_node),err,error,*999)
    ENDDO
    CALL nodes_create_finish(nodes,err,error,*999)
    CALL checked_deallocate( user_nodal_number_map_global_nodal_number )

    !IF(ALLOCATED(USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER)) DEALLOCATE(USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER)
    !ALLOCATE(USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER(NUMBER_OF_NODES), STAT=ERR)
    !IF(ERR/=0) CALL FlagError("can not allocate nodal number mapping for output",ERR,ERROR,*999)
    !USER_NODAL_NUMBER_MAP_GLOBAL_NODAL_NUMBER(:)=LIST_NODAL_NUMBER(:)
    !IF(ALLOCATED(LIST_NODAL_NUMBER)) DEALLOCATE(LIST_NODAL_NUMBER)

    !calculate the number of mesh components
    IF(master_computational_number==my_computational_node_number) THEN

      !MESH_COMPONENT_LOOKUP is used to store the difference between field components in term of basis property.
      CALL reallocate_2d( mesh_component_lookup, number_of_components, number_of_components, &
        & "can not allocate list of mesh components", err, error, *999 )

      CALL reallocate( list_str, number_of_components, &
        & "can not allocate list of str", err, error, *999 )
      !initialize MESH_COMPONENT_LOOKUP and assume each field component has the same mesh component

      DO idx_comp=1,number_of_components
        mesh_component_lookup(idx_comp,idx_comp)=1
      ENDDO
      idx_exelem=0

      !checking field component's mesh component by checking the basis
      DO WHILE(idx_exelem<number_of_exelem_files)

        file_id=1030+idx_exelem
        !checking the next file
        file_name=name//".part"//trim(number_to_vstring(idx_exelem,"*",err,error))//".exelem"
        CALL field_io_fortran_file_open(file_id, file_name, file_status, err,error,*999)
        field_section=.false.
        file_end=.false.

        DO WHILE(.NOT.file_end)
          CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)

          !check the beginning of element section
          IF((.NOT.start_of_element_section).AND.(verify(cmiss_keyword_shape,line)==0)) THEN
            start_of_element_section=.true.
          ENDIF

          !check whether it is a new header for another group of elements
          IF(start_of_element_section.AND.(verify(cmiss_keyword_fields,line)==0)) THEN

            !collect header information
            pos=index(line,cmiss_keyword_fields)
            line=remove(line,1, pos+len_trim(cmiss_keyword_fields)-1)
            number_of_fields=string_to_integer(line, err,error)
            idx_comp=0
            DO idx_field=1,number_of_fields
              CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
              pos=index(line,cmiss_keyword_components)
              line=remove(line,1, pos+len_trim(cmiss_keyword_components)-1)
              number_of_comp=string_to_integer(line, err,error)
              DO idx_comp1=1,number_of_comp
                idx_comp=idx_comp+1
                CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
                pos=index(line,".")
                line=remove(line,1, pos)
                line=trim(adjustl(line))
                pos=index(line,",")
                line=remove(line,pos, len(line))
                list_str(idx_comp)= line
                CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
                pos=index(line, cmiss_keyword_nodes)
                line=remove(line,1, pos+len_trim(cmiss_keyword_nodes)-1)
                number_of_node=string_to_integer(line, err,error)
                DO idx_node1=1, number_of_node*number_nodal_lines
                  CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
                ENDDO !idx_node1
              ENDDO !idx_comp1
            ENDDO !idx_field

            !compare the bases. Since the geometrical topology is same, if the bases are the same in the same topology,
            !they should be in the same mesh component
            IF(sum(mesh_component_lookup)<number_of_components*number_of_components) THEN
              DO idx_comp=1, number_of_components
                DO idx_comp1=idx_comp+1, number_of_components
                  IF(mesh_component_lookup(idx_comp1,idx_comp)==0) THEN
                    IF(list_str(idx_comp1)/=list_str(idx_comp)) THEN
                      mesh_component_lookup(idx_comp1,idx_comp)=1
                      mesh_component_lookup(idx_comp,idx_comp1)=1
                    ENDIF
                  ENDIF
                ENDDO !idx_comp1
              ENDDO !idx_comp
            ELSE
              idx_exelem=number_of_exelem_files !jump out of the loop
            ENDIF
          ENDIF !START_OF_ELEMENT_SECTION==.TRUE..AND.VERIFY(CMISS_KEYWORD_FIELDS,LINE)==0

          !!check whether they have same numbers of field components
          !IF(FIELD_SECTION==.TRUE..AND.START_OF_ELEMENT_SECTION==.TRUE..AND.(VERIFY(CMISS_KEYWORD_COMPONENTS,LINE)==0)) THEN
          !  IF(idx_field>=NUMBER_OF_FIELDS) THEN
          !    FIELD_SECTION=.FALSE.
          !  ENDIF
          !ENDIF !FIELD_SECTION==.TRUE..AND.START_OF_ELEMENT_SECTION==.TRUE..AND.VERIFY(CMISS_KEYWORD_COMPONENTS,LINE
        ENDDO !(FILE_END==.FALSE.)
        CALL field_io_fortran_file_close(file_id, err,error,*999)
        idx_exelem=idx_exelem+1
      ENDDO !idx_exelem=0, NUMBER_OF_EXELEM_FILES-1

      !calculate the number of mesh components
      CALL reallocate( mesh_components_of_field_components, number_of_components, &
          & "can not allocate list of field components", err, error, *999 )

      DO idx_comp=1, number_of_components
        mesh_components_of_field_components(idx_comp)=idx_comp
      ENDDO
      DO idx_comp=1, number_of_components
        DO idx_comp1=idx_comp+1, number_of_components
          IF(mesh_components_of_field_components(idx_comp)==idx_comp) THEN
            IF(mesh_component_lookup(idx_comp1,idx_comp)==0) mesh_components_of_field_components(idx_comp1)=idx_comp
          ENDIF
        ENDDO !idx_comp1
      ENDDO !idx_comp
      CALL checked_deallocate( mesh_component_lookup )
      number_of_mesh_components=0
      idx_comp1=0
      DO idx_comp=1,number_of_components
        IF(mesh_components_of_field_components(idx_comp)==idx_comp) THEN
          idx_comp1=idx_comp1+1
          mesh_components_of_field_components(idx_comp)=idx_comp1
          number_of_mesh_components=number_of_mesh_components+1
        ENDIF
      ENDDO
      !IF(ALLOCATED(MESH_COMPONENTS_OF_FIELD_COMPONENTS)) DEALLOCATE(MESH_COMPONENTS_OF_FIELD_COMPONENTS)
      !ALLOCATE(MESH_COMPONENTS_OF_FIELD_COMPONENTS(NUMBER_OF_COMPONENTS),STAT=ERR)
      !IF(ERR/=0) CALL FlagError("can not allocate mesh components of field components for output",ERR,ERROR,*999)
      !MESH_COMPONENTS_OF_FIELD_COMPONENTS(:)=LIST_FIELD_COMPONENTS(:)
      !ALLOCATE(LIST_MESH_COMPONENTS(NUMBER_OF_MESH_COMPONENTS),STAT=ERR)
      !IF(ERR/=0) CALL FlagError("ALLOCATEDcan not allocate list of mesh components",ERR,ERROR,*999)
      !idx_comp1=0
      !DO idx_comp=1,NUMBER_OF_COMPONENTS
      !   IF(LIST_FIELD_COMPONENTS(idx_comp)==idx_comp) THEN
      !        idx_comp1=idx_comp1+1
      !      LIST_MESH_COMPONENTS(idx_comp1)=idx_comp
      !   ENDIF
      !ENDDO
    ENDIF !MASTER_COMPUTATIONAL_NUMBER==my_computational_node_number

    !broadcasting the number of mesh components
    CALL mpi_bcast(number_of_components,1,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
    CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
    CALL mpi_bcast(number_of_mesh_components,1,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
    CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
    CALL mesh_number_of_components_set(mesh,number_of_mesh_components,err,error,*999)

    CALL reallocate( elements_ptr, number_of_mesh_components, &
      & "can not allocate list of mesh element pointers", err, error, *999 )

    IF(basis_functions%NUMBER_BASIS_FUNCTIONS<=0)  THEN
      CALL basis_create_start(1,basis,err,error,*999)
      CALL basis_number_of_xi_set(basis,number_of_dimensions,err,error,*999)
      CALL basis_create_finish(basis,err,error,*999)
    ENDIF

    DO idx_comp=1, number_of_mesh_components
      CALL mesh_topology_elements_create_start(mesh,idx_comp,basis_functions%BASES(1)%PTR,elements_ptr(idx_comp)%PTR, &
          & err,error,*999)
    ENDDO

    !Collect the elemental numberings (elemental labels)
    CALL reallocate( list_element_number, number_of_elements, &
        & "can not allocate list of elemental number", err, error, *999 )

    IF(master_computational_number==my_computational_node_number) THEN
      !the file name has to start from zero in a ascended order without break
      idx_elem=1
      DO idx_exelem=0, number_of_exelem_files-1

         file_id=1030+idx_exelem
         file_name=name//".part"//trim(number_to_vstring(idx_exelem,"*",err,error))//".exelem"
         CALL field_io_fortran_file_open(file_id, file_name, file_status, err,error,*999)
         start_of_element_section=.false.
         file_end=.false.


         DO WHILE(.NOT.file_end)
           CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)

           IF((.NOT.start_of_element_section).AND.(verify(cmiss_keyword_shape,line)==0)) THEN
             start_of_element_section=.true.
           ENDIF

           IF(start_of_element_section.AND.(verify(cmiss_keyword_element,line)==0)) THEN
             pos=index(line,cmiss_keyword_element)
             line=remove(line,1, pos+len_trim(cmiss_keyword_element)-1)
             shape_index(:)=0
             CALL string_to_muti_integers_vs(line, shape_size, shape_index(:), err,error, *999)
             IF(number_of_dimensions==3) THEN
               list_element_number(idx_elem)=shape_index(1)
             ELSE IF(number_of_dimensions==2) THEN
               list_element_number(idx_elem)=shape_index(2)
             ELSE IF(number_of_dimensions==1) THEN
               list_element_number(idx_elem)=shape_index(3)
             ELSE
               CALL flagerror("Non recognized dimension size during reading elemental numbering",err,error,*999)
             ENDIF
             idx_elem=idx_elem+1
           ENDIF
         ENDDO !(FILE_END==.FALSE.)

         CALL field_io_fortran_file_close(file_id, err,error,*999)
      ENDDO !idx_exelem=0
      CALL list_sort(list_element_number, err, error, *999)
    ENDIF !MASTER_COMPUTATIONAL_NUMBER==my_computational_node_number

    !broadcast the list of elements for mapping gloabl numbers and user numbers (elemental labels)
    CALL mpi_bcast(list_element_number,number_of_elements,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
    CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
    !change the mapping between global elemental numbering and user elemental numbering

    DO idx_elem=1,number_of_elements
      DO idx_comp=1, number_of_mesh_components
        IF(idx_elem/=list_element_number(idx_elem)) &
            & CALL meshelements_elementusernumberset(idx_elem,list_element_number(idx_elem), &
            & elements_ptr(idx_comp)%PTR,err,error,*999)
      ENDDO
    ENDDO

    !creating topological information for each mesh component
    CALL mpi_bcast(number_of_exelem_files,1,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
    CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
    !ALLOCATE(LIST_BASES(NUMBER_OF_COMPONENTS),STAT=ERR)
    !IF(ERR/=0) CALL FlagError("can not allocate list of bases",ERR,ERROR,*999)

    CALL reallocate( list_comp_nodes, number_of_components, &
        & "Could not allocate list of component nodal index ", err, error, *999 )

    file_end=.true.
    idx_exelem=-1
    file_id=1030

    DO WHILE(idx_exelem<number_of_exelem_files)

      CALL mpi_bcast(file_end,1,mpi_logical,master_computational_number,mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)

      IF(file_end) THEN
        idx_exelem=idx_exelem+1
        INQUIRE(unit=file_id, opened=file_open)
        IF(file_open) CALL field_io_fortran_file_close(file_id, err,error,*999)
        IF(idx_exelem>=number_of_exelem_files) EXIT
      ENDIF

      !goto the start of mesh part
      IF(master_computational_number==my_computational_node_number) THEN

        !IF(FILE_END==.FALSE..AND.START_OF_ELEMENT_SECTION==.FALSE.) THEN
        !  !check the beginning of element section
        !  DO WHILE(VERIFY(CMISS_KEYWORD_SHAPE,LINE)/=0)
        !    CALL FIELD_IO_FORTRAN_FILE_READ_STRING(FILE_ID, LINE, FILE_END, ERR,ERROR, *999)
        !  ENDDO
        !  START_OF_ELEMENT_SECTION=.TRUE.
        !ENDIF

        IF(file_end) THEN
          file_id=1030+idx_exelem
          !checking the next file
          file_name=name//".part"//trim(number_to_vstring(idx_exelem,"*",err,error))//".exelem"
          CALL field_io_fortran_file_open(file_id, file_name, file_status, err,error,*999)
          file_end=.false.
          section_start=.false.
          start_of_element_section=.false.
        ENDIF


        IF((.NOT.file_end).AND.(.NOT.start_of_element_section)) THEN !..AND.SECTION_START==.FALSE.) THEN
          !find a new header
          DO WHILE(verify(cmiss_keyword_scale_factor_sets,line)/=0)
            CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
          ENDDO
          start_of_element_section=.true.
        ENDIF

        !have not touched the end
        IF((.NOT.file_end).AND.start_of_element_section.AND.(.NOT.section_start)) THEN
          section_start=.true.
          !CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,LINE,ERR,ERROR,*999)
          !CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,FILE_NAME,ERR,ERROR,*999)

          pos=index(line,cmiss_keyword_scale_factor_sets)
          line=remove(line,1, pos+len_trim(cmiss_keyword_scale_factor_sets)-1)
          number_of_scalesets=string_to_integer(line, err,error)
          idx_mesh_comp=1
          !skip factors
          number_scaling_factor_lines=0
          DO idx_scl=1,number_of_scalesets
            CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
            pos=index(line,cmiss_keyword_scale_factors)
            line=remove(line,1, pos+len_trim(cmiss_keyword_scale_factors)-1)
            num_scl=string_to_integer(line, err,error)
            num_scl_line=num_scl/number_scaling_factors_in_line
            IF(num_scl_line*number_scaling_factors_in_line/=num_scl) num_scl_line=num_scl_line+1
            number_scaling_factor_lines=number_scaling_factor_lines+num_scl_line
          ENDDO
          !skip nodes line
          CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
          pos=index(line,cmiss_keyword_nodes)
          line=remove(line,1, pos+len_trim(cmiss_keyword_nodes)-1)
          number_of_node=string_to_integer(line, err,error)

          CALL reallocate( list_elemental_nodes, number_of_node, &
            & "Could not allocate list of elemental nodes", err, error, *999 )
          !CALL WRITE_STRING_VALUE(GENERAL_OUTPUT_TYPE,"number_of_node:",number_of_node,ERR,ERROR,*999)

          CALL reallocate_2d( list_comp_nodal_index, number_of_components,number_of_node, &
            & "Could not allocate list of component nodal index ", err, error, *999 )

          !read the header
          CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
          pos=index(line,cmiss_keyword_fields)
          line=remove(line,1, pos+len_trim(cmiss_keyword_fields)-1)
          number_of_fields=string_to_integer(line, err,error)
          idx_comp=0
          DO idx_field=1,number_of_fields
            CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
            pos=index(line,cmiss_keyword_components)
            line=remove(line,1, pos+len_trim(cmiss_keyword_components)-1)
            number_of_comp=string_to_integer(line, err,error)
            DO idx_comp1=1,number_of_comp
              idx_comp=idx_comp+1
              CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
              pos=index(line,".")
              line=remove(line,1, pos)
              line=trim(adjustl(line))
              pos=index(line,",")
              line=remove(line,pos, len(line))
              list_str(idx_comp)= line
              CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
              pos=index(line, cmiss_keyword_nodes)
              line=remove(line,1, pos+len_trim(cmiss_keyword_nodes)-1)
              number_of_node=string_to_integer(line, err,error)
              list_comp_nodes(idx_comp)=number_of_node
              DO idx_node1=1, number_of_node
                CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
                pos=index(line, ".")
                line=remove(line,pos, len(line))
                list_comp_nodal_index(idx_comp,idx_node1)=string_to_integer(line, err,error)
                CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
                CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
              ENDDO !idx_node1
            ENDDO !idx_comp1
          ENDDO !idx_field

          CALL reallocate_2d( interpolation_xi, number_of_components,number_of_dimensions, &
            & "Could not allocate list of interpolation types", err, error, *999 )

          CALL field_io_fill_basis_info(interpolation_xi, list_str, number_of_components, err,error,*999)
          CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
          !CALL WRITE_STRING(GENERAL_OUTPUT_TYPE,LINE,ERR,ERROR,*999)
        ENDIF  !FILE_END==.FALSE..AND.START_OF_ELEMENT_SECTION==.TRUE..AND.SECTION_START==.TRUE.

        IF((.NOT.file_end).AND.start_of_element_section.AND.section_start) THEN

          pos=index(line,cmiss_keyword_element)
          line=remove(line,1, pos+len_trim(cmiss_keyword_element)-1)
          CALL string_to_muti_integers_vs(line, shape_size, shape_index(:), err,error, *999)


          DO WHILE(verify(cmiss_keyword_node,line)/=0)
            CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
          ENDDO

          CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
          CALL string_to_muti_integers_vs(line, number_of_node, list_elemental_nodes, err, error, *999)

          !skip scaling factors
          CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
          DO idx_scl=1, number_scaling_factor_lines
            CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
          ENDDO

          IF(.NOT.file_end) THEN
            CALL field_io_fortran_file_read_string(file_id, line, file_end, err,error, *999)
            IF(verify(cmiss_keyword_scale_factor_sets,line)==0) section_start=.true.
          ENDIF
          !CALL WRITE_STRING_VALUE(GENERAL_OUTPUT_TYPE,"FILE_END:",FILE_END,ERR,ERROR,*999)
        ENDIF
      ENDIF !MASTER_COMPUTATIONAL_NUMBER

      CALL mpi_bcast(number_of_node,1,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
      !CALL WRITE_STRING_VALUE(GENERAL_OUTPUT_TYPE,"SIZE LIST_ELEMENTAL_NODES:",SIZE(LIST_ELEMENTAL_NODES),ERR,ERROR,*999)

      IF(master_computational_number/=my_computational_node_number) THEN
        CALL reallocate( list_elemental_nodes, number_of_node, &
          & "Could not allocate list of elemental nodes", err, error, *999 )
        CALL reallocate_2d( list_comp_nodal_index, number_of_components, number_of_node, &
          & "Could not allocate list of component nodal index ", err, error, *999 )
        CALL reallocate_2d( interpolation_xi, number_of_components,number_of_dimensions, &
          & "Could not allocate list of interpolation types", err, error, *999 )
        CALL reallocate( mesh_components_of_field_components, number_of_components, &
          & "Could not allocate list of mesh components of field", err, error, *999 )
        !DO idx_comp=1,NUMBER_OF_COMPONENTS
        !  DO idx_dim=1,NUMBER_OF_DIMENSIONS
        !    IF(idx_dim==NUMBER_OF_DIMENSIONS) THEN
        !      LINE=LIST_STR(idx_comp)
        !    ELSE
        !      pos=INDEX(LINE,"*")
        !      LINE=EXTRACT(LIST_STR(idx_comp),1, pos-1)
        !      LIST_STR(idx_comp)=REMOVE(LIST_STR(idx_comp),1,pos)
        !    ENDIF
        !    CALL FieldIO_TranslateLabelIntoInterpolationType(INTERPOLATION_XI(idx_comp, idx_dim), LINE, ERR, ERROR, *999)
        !  ENDDO
        !ENDDO
      ENDIF !MASTER_COMPUTATIONAL_NUMBER/=my_computational_node_number

      !CALL WRITE_STRING_VALUE(GENERAL_OUTPUT_TYPE,"LIST_ELEMENTAL_NODES:",LIST_ELEMENTAL_NODES(1),ERR,ERROR,*999)
      CALL mpi_bcast(list_elemental_nodes,number_of_node,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
      CALL mpi_bcast(list_comp_nodal_index,number_of_node*number_of_components,mpi_integer,master_computational_number, &
          & mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
      CALL mpi_bcast(shape_index,shape_size,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
      CALL mpi_bcast(list_comp_nodes,number_of_components,mpi_integer,master_computational_number,mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
      CALL mpi_bcast(mesh_components_of_field_components,number_of_components,mpi_integer,master_computational_number, &
          & mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
      CALL mpi_bcast(interpolation_xi,number_of_components*number_of_dimensions,mpi_integer,master_computational_number,&
          & mpi_comm_world,mpi_ierror)
      CALL mpi_error_check("MPI_BCAST",mpi_ierror,err,error,*999)
      !CALL WRITE_STRING_VALUE(GENERAL_OUTPUT_TYPE,"LIST_ELEMENTAL_NODES:",LIST_ELEMENTAL_NODES(1),ERR,ERROR,*999)
      current_mesh_comp=1
      DO idx_comp=1, number_of_components
        !CALL WRITE_STRING_VALUE(GENERAL_OUTPUT_TYPE,"LIST_ELEMENTAL_NODES:",LIST_ELEMENTAL_NODES(1),ERR,ERROR,*999)
        IF(number_of_dimensions==3) THEN
          CALL list_search(list_element_number, shape_index(1),global_element_number, err,error,*999)
        ELSE IF(number_of_dimensions==2) THEN
          CALL list_search(list_element_number, shape_index(2),global_element_number, err,error,*999)
        ELSE IF(number_of_dimensions==1) THEN
          CALL list_search(list_element_number, shape_index(3),global_element_number, err,error,*999)
        ELSE
          CALL flagerror("Non recognized dimension size during reading elemental numbering",err,error,*999)
        ENDIF

        IF(mesh_components_of_field_components(idx_comp)==current_mesh_comp) THEN
          !find out whether the basis has been created
          pos=0
          DO idx_basis=1, basis_functions%NUMBER_BASIS_FUNCTIONS
            IF(sum(basis_functions%BASES(idx_basis)%PTR%INTERPOLATION_XI(:)-interpolation_xi(idx_comp,:))==0) THEN
              pos=idx_basis
              EXIT
            ENDIF
          ENDDO

          IF(pos==0) THEN
            IF(ASSOCIATED(basis)) NULLIFY(basis)
            CALL basis_create_start(basis_functions%NUMBER_BASIS_FUNCTIONS+1,basis,err,error,*999)
            CALL basis_number_of_xi_set(basis,number_of_dimensions,err,error,*999)
            CALL basis_interpolation_xi_set(basis,interpolation_xi(idx_comp,:),err,error,*999)
            CALL basis_create_finish(basis,err,error,*999)
          ENDIF

          CALL mesh_topology_elements_element_basis_set(global_element_number,elements_ptr( &
            & mesh_components_of_field_components(idx_comp))%PTR,basis,err,error,*999)
          !CALL WRITE_STRING_VALUE(GENERAL_OUTPUT_TYPE,"LIST_ELEMENTAL_NODES:",LIST_ELEMENTAL_NODES(1),ERR,ERROR,*999)
          CALL mesh_topology_elements_element_nodes_set(global_element_number,elements_ptr( &
            & mesh_components_of_field_components(idx_comp))%PTR,list_elemental_nodes(list_comp_nodal_index(idx_comp,:)), &
            & err,error,*999)
          !CALL WRITE_STRING_VALUE(GENERAL_OUTPUT_TYPE,"LIST_ELEMENTAL_NODES:",LIST_ELEMENTAL_NODES(1),ERR,ERROR,*999)
          current_mesh_comp=current_mesh_comp+1
        ENDIF
      ENDDO !idx_comp
      !CALL WRITE_STRING_VALUE(GENERAL_OUTPUT_TYPE,"FILE_END:",FILE_END,ERR,ERROR,*999)

    ENDDO !idx_exelem<NUMBER_OF_EXELEM_FILES

    DO idx_comp=1, number_of_mesh_components
      CALL mesh_topology_elements_create_finish(elements_ptr(idx_comp)%PTR, err,error,*999)
    ENDDO
    CALL mesh_create_finish(mesh,err,error,*999)

    CALL checked_deallocate( list_element_number )
    CALL checked_deallocate( elements_ptr )
    !IF(ALLOCATED(LIST_NODAL_NUMBER)) DEALLOCATE(LIST_NODAL_NUMBER)
    CALL checked_deallocate( list_elemental_nodes )
    CALL checked_deallocate( list_str )
    CALL checked_deallocate( interpolation_xi )
    !IF(ALLOCATED(LIST_FIELD_COMPONENTS)) DEALLOCATE(LIST_FIELD_COMPONENTS)
    CALL checked_deallocate( mesh_component_lookup )
    CALL checked_deallocate( list_comp_nodal_index )
    CALL checked_deallocate( list_comp_nodes )
    CALL checked_deallocate( user_nodal_number_map_global_nodal_number )
    IF(ASSOCIATED(basis)) NULLIFY(basis)

    exits("FIELD_IO_IMPORT_GLOBAL_MESH")
    RETURN
999 errorsexits("FIELD_IO_IMPORT_GLOBAL_MESH",err,error)
    RETURN 1
  END SUBROUTINE field_io_import_global_mesh

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

  SUBROUTINE fieldio_translatelabelintointerpolationtype(INTERPOLATION, LABEL_TYPE, ERR, ERROR, *)
    !Argument variables
    INTEGER(INTG), INTENT(INOUT) :: INTERPOLATION
    TYPE(varying_string), INTENT(IN) :: LABEL_TYPE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    SELECT CASE(char(label_type))
    CASE("l.Lagrange")
      interpolation=basis_linear_lagrange_interpolation
    CASE("q.Lagrange")
      interpolation=basis_quadratic_lagrange_interpolation
    CASE("c.Lagrange")
      interpolation=basis_cubic_lagrange_interpolation
    CASE("c.Hermite")
      interpolation=basis_cubic_hermite_interpolation
    CASE("LagrangeHermite")
      interpolation=basis_quadratic1_hermite_interpolation
    CASE("HermiteLagrange")
      interpolation=basis_quadratic2_hermite_interpolation
    CASE DEFAULT
      CALL flagerror("Invalid interpolation type",err,error,*999)
    END SELECT

    exits("FieldIO_TranslateLabelIntoInterpolationType")
    RETURN
999 errorsexits("FieldIO_TranslateLabelIntoInterpolationType",err,error)
    RETURN 1
    
  END SUBROUTINE fieldio_translatelabelintointerpolationtype

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

  SUBROUTINE fieldio_calculatesimplexscaleandnodecounts(BASIS, num_scl, num_node, ERR, ERROR, * )
    !Argument variables
    TYPE(basis_type), INTENT(IN) :: BASIS
    INTEGER(INTG), INTENT(INOUT) :: num_scl
    INTEGER(INTG), INTENT(INOUT) :: num_node
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    !Local variables
    INTEGER(INTG) :: n

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

    IF(basis%NUMBER_OF_XI==0) CALL flagerror("number of xi in the basis is zero",err,error,*999)

    n = basis%NUMBER_OF_XI

    !Simplex-type interpolations must be the same in all xi.
    SELECT CASE(basis%INTERPOLATION_XI(1))
      CASE(basis_linear_simplex_interpolation)
        num_node = n + 1
      CASE(basis_quadratic_simplex_interpolation)
        num_node = ( n + 1 ) * ( n + 2 ) / 2
      CASE(basis_cubic_simplex_interpolation)
        num_node = ( n + 1 ) * ( n + 2 ) * ( n + 3 ) / 6
      CASE DEFAULT
        CALL flagerror( "Invalid interpolation type", err, error, *999 )
    END SELECT

    num_scl = num_node

    exits("FieldIO_CalculateSimplexScaleAndNodeCounts")
    RETURN
999 errorsexits("FieldIO_CalculateSimplexScaleAndNodeCounts",err,error)
    RETURN 1
  END SUBROUTINE fieldio_calculatesimplexscaleandnodecounts

  !
  !================================================================================================================================
  !
  SUBROUTINE field_io_calculate_tp_scale_and_node_counts(BASIS, num_scl, num_node, ERR, ERROR, * )
    !Argument variables
    TYPE(basis_type), INTENT(IN) :: BASIS
    INTEGER(INTG), INTENT(INOUT) :: num_scl
    INTEGER(INTG), INTENT(INOUT) :: num_node
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: ni

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

    IF(basis%NUMBER_OF_XI==0) CALL flagerror("number of xi in the basis is zero",err,error,*999)

    num_scl=1;
    num_node=1
    DO ni=1,basis%NUMBER_OF_XI
      SELECT CASE(basis%INTERPOLATION_XI(ni))
      CASE(basis_linear_lagrange_interpolation)
        num_scl=num_scl*2
        num_node=num_node*2
      CASE(basis_quadratic_lagrange_interpolation)
        num_scl=num_scl*3
        num_node=num_node*3
      CASE(basis_cubic_lagrange_interpolation)
        num_scl=num_scl*4
        num_node=num_node*4
      CASE(basis_cubic_hermite_interpolation)
        num_scl=num_scl*4
        num_node=num_node*2
      CASE(basis_quadratic1_hermite_interpolation)
        num_scl=num_scl*3
        num_node=num_node*2
      CASE(basis_quadratic2_hermite_interpolation)
        num_scl=num_scl*3
        num_node=num_node*2
      CASE DEFAULT
        CALL flagerror( "Invalid interpolation type", err, error, *999 )
      END SELECT
    ENDDO !ni

    exits("FIELD_IO_CALCULATE_TP_SCALE_AND_NODE_COUNTS")
    RETURN
999 errorsexits("FIELD_IO_CALCULATE_TP_SCALE_AND_NODE_COUNTS",err,error)
    RETURN 1
  END SUBROUTINE field_io_calculate_tp_scale_and_node_counts

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

  FUNCTION findmylocaldomainnumber( mapping, myComputationalNodeNumber )
    TYPE(domain_global_mapping_type) :: mapping
    INTEGER(INTG), INTENT(IN) :: myComputationalNodeNumber

    INTEGER(INTG) :: FindMyLocalDomainNumber

    INTEGER(INTG) :: domainIndex
    INTEGER(INTG) :: myDomainIndex

    DO domainindex = 1, mapping%NUMBER_OF_DOMAINS
      IF( mapping%DOMAIN_NUMBER( domainindex ) == mycomputationalnodenumber ) THEN
        mydomainindex = domainindex
        EXIT
      ENDIF
    ENDDO

    findmylocaldomainnumber = mapping%LOCAL_NUMBER( mydomainindex )
  END FUNCTION findmylocaldomainnumber

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

  SUBROUTINE fieldio_exportelementalgroupheaderfortran( global_number, MAX_NODE_COMP_INDEX,NUM_OF_SCALING_FACTOR_SETS, &
    & list_comp_scale, my_computational_node_number, elementalinfoset, sessionhandle, err,error, *)
    !Argument variables
    INTEGER(INTG), INTENT(IN) :: global_number
    INTEGER(INTG), INTENT(INOUT) ::  MAX_NODE_COMP_INDEX
    INTEGER(INTG), INTENT(INOUT) :: NUM_OF_SCALING_FACTOR_SETS
    INTEGER(INTG), INTENT(INOUT) :: LIST_COMP_SCALE(:)
    INTEGER(INTG), INTENT(IN) :: my_computational_node_number
    TYPE(field_io_component_info_set), INTENT(INOUT) :: elementalInfoSet
    INTEGER(INTG), INTENT(IN) :: sessionHandle
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: i,LENGTH
    INTEGER(INTG) :: NUMBER_OF_UNIQUE_NODES
    CHARACTER(LEN=MAXSTRLEN) :: fvar_name
    CHARACTER(LEN=1, KIND=C_CHAR) :: cvar_name(maxstrlen+1)
    TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
    TYPE(field_ptr_type), ALLOCATABLE :: listScaleFields(:)
    TYPE(field_variable_type), POINTER :: variable_ptr
    TYPE(domain_type), POINTER :: componentDomain !The domain mapping to calculate nodal mappings
    TYPE(domain_elements_type), POINTER :: DOMAIN_ELEMENTS ! domain nodes
    TYPE(domain_element_type), POINTER :: MAX_NODE_ELEMENT
    TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES,MAX_ELEMENT_DOMAIN_NODES ! domain nodes
    TYPE(basis_type), POINTER :: BASIS
    TYPE(basis_ptr_type), ALLOCATABLE :: listScaleBases(:)
    TYPE(field_variable_component_type), POINTER :: component
    INTEGER(INTG), ALLOCATABLE :: GROUP_LOCAL_NUMBER(:), GROUP_SCALE_FACTORS(:)
    INTEGER(INTG), ALLOCATABLE :: GROUP_NODE(:), GROUP_VARIABLES(:)
    !INTEGER(C_INT), TARGET :: INTERPOLATION_XI(3),ELEMENT_DERIVATIVES(64*64),NUMBER_OF_DERIVATIVES(64), NODE_INDEXES(128)
    INTEGER(C_INT), ALLOCATABLE, TARGET :: INTERPOLATION_XI(:),ELEMENT_DERIVATIVES(:),NUMBER_OF_DERIVATIVES(:), NODE_INDEXES(:)
    INTEGER(C_INT), ALLOCATABLE, TARGET :: SCALE_INDEXES(:) !Array for holding scale indexes, useful for collapsed nodes.
    INTEGER(INTG) :: nn, nx, ny, nz, NodesX, NodesY, NodesZ, mm, NUM_OF_VARIABLES, MAX_NUM_NODES !NUM_OF_NODES
    INTEGER(INTG) :: local_number, interpType, NODE_NUMBER, NODE_NUMBER_COUNTER, NODE_NUMBER_COLLAPSED, NUMBER_OF_ELEMENT_NODES
    INTEGER(INTG) :: num_scl, num_node, comp_idx, scaleIndex, scaleIndex1, var_idx, derivativeIndex !value_idx field_idx global_var_idx comp_idx1 ny2
    INTEGER(INTG) :: NODE_LOCAL_NUMBER,NODE_USER_NUMBER,MAX_ELEMENT_LOCAL_NUMBER,MAX_ELEMENT_USER_NUMBER
    LOGICAL :: SAME_SCALING_SET

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

    !SANDER
    ALLOCATE(interpolation_xi(3), stat = err)

    !colllect nodal header information for IO first

    !collect maximum number of nodal derivatives, number of fields and variables
    num_of_scaling_factor_sets=0
    num_of_variables=0
    max_num_nodes=0
    max_node_comp_index=0
    NULLIFY(variable_ptr)

    CALL reallocate( group_local_number, elementalinfoset%NUMBER_OF_COMPONENTS, &
        & "Could not allocate GROUP_LOCAL_NUMBER in exelem header", err, error, *999 )
    CALL reallocate( listscalebases, elementalinfoset%NUMBER_OF_COMPONENTS, &
        & "Could not allocate listScaleBases in exelem header", err, error, *999 )
    CALL reallocate( listscalefields, elementalinfoset%NUMBER_OF_COMPONENTS, &
        & "Could not allocate listScaleFields in exelem header", err, error, *999 )

    !collect scale factor information
    DO comp_idx=1,elementalinfoset%NUMBER_OF_COMPONENTS
      !calculate the number of variables
      IF (.NOT.ASSOCIATED(variable_ptr, TARGET=elementalinfoset%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE)) THEN
        num_of_variables=num_of_variables+1
        variable_ptr=>elementalinfoset%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE
      ENDIF

      !finding the local numbering through the global to local mapping
      componentdomain=>elementalinfoset%COMPONENTS(comp_idx)%PTR%DOMAIN
      !get the domain index for this variable component according to my own computional node number
      local_number = findmylocaldomainnumber( componentdomain%MAPPINGS%ELEMENTS%GLOBAL_TO_LOCAL_MAP( global_number ),&
          & my_computational_node_number )
      group_local_number(comp_idx)=local_number
      !use local domain information find the out the maximum number of derivatives
      domain_elements=>componentdomain%TOPOLOGY%ELEMENTS
      domain_nodes=>componentdomain%TOPOLOGY%NODES
      basis=>domain_elements%ELEMENTS(local_number)%BASIS
      IF(basis%NUMBER_OF_NODES>max_num_nodes) THEN
        max_node_comp_index=comp_idx
        max_node_element => domain_elements%ELEMENTS(local_number)
        max_num_nodes=basis%NUMBER_OF_NODES
        max_element_domain_nodes=>componentdomain%TOPOLOGY%NODES
      ENDIF
      !IF(.NOT.BASIS%DEGENERATE)  THEN
      IF(comp_idx == 1) THEN
        num_of_scaling_factor_sets = num_of_scaling_factor_sets + 1
        listscalebases( num_of_scaling_factor_sets )%PTR => basis
        listscalefields( num_of_scaling_factor_sets )%PTR => variable_ptr%FIELD
        list_comp_scale(comp_idx)=num_of_scaling_factor_sets
      ELSE
        same_scaling_set=.false.
        DO scaleindex1=1, num_of_scaling_factor_sets
          IF(basis%GLOBAL_NUMBER == listscalebases(scaleindex1)%PTR%GLOBAL_NUMBER) THEN
            IF(variable_ptr%FIELD%SCALINGS%SCALING_TYPE /= listscalefields(scaleindex1)%PTR%SCALINGS%SCALING_TYPE) THEN
              CALL flag_warning("Fields "//trim(number_to_vstring(listscalefields(scaleindex1)%PTR%USER_NUMBER,"*",err,error))// &
                  & " and "//trim(number_to_vstring(variable_ptr%FIELD%USER_NUMBER,"*",err,error))// &
                  & " have components that use basis number "//trim(number_to_vstring(basis%GLOBAL_NUMBER,"*",err,error))// &
                  & " but have different scaling types. ",err,error,*999)
            ENDIF
            IF(variable_ptr%FIELD%SCALINGS%SCALING_TYPE == listscalefields(scaleindex1)%PTR%SCALINGS%SCALING_TYPE) THEN
              same_scaling_set=.true.
              list_comp_scale(comp_idx)=scaleindex1
              EXIT
            ENDIF
          ENDIF
        ENDDO !scaleIndex1
        IF(.NOT.same_scaling_set) THEN
          num_of_scaling_factor_sets=num_of_scaling_factor_sets+1
          listscalebases( num_of_scaling_factor_sets )%PTR => basis
          listscalefields( num_of_scaling_factor_sets )%PTR => variable_ptr%FIELD
          list_comp_scale(comp_idx)=num_of_scaling_factor_sets
        ENDIF
      ENDIF
      !ENDIF !BASIS%DEGENERATE=.FALSE.
    ENDDO !comp_idx
    !!Allocate the memory for group of field components
    CALL reallocate( group_variables, num_of_variables, &
        & "Could not allocate temporary variable buffer in IO", err, error, *999 )

    !!Allocate the memory for group of maximum number of derivatives
    CALL reallocate( group_scale_factors, num_of_scaling_factor_sets, &
        & "Could not allocate temporary variable buffer in IO", err, error, *999 )

    CALL reallocate( group_node, num_of_scaling_factor_sets, &
        & "Could not allocate temporary variable buffer in IO", err, error, *999 )

    !fill information into the group of fields and variables
    NULLIFY(variable_ptr)
    num_of_variables=0
    DO comp_idx=1,elementalinfoset%NUMBER_OF_COMPONENTS
      !calculate the number of variables
      IF (.NOT.ASSOCIATED(variable_ptr, TARGET=elementalinfoset%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE)) THEN
        num_of_variables=num_of_variables+1
        variable_ptr=>elementalinfoset%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE
      ENDIF
      group_variables(num_of_variables)=group_variables(num_of_variables)+1
    ENDDO  !comp_idx

    DO scaleindex = 1, num_of_scaling_factor_sets
      basis => listscalebases( scaleindex )%PTR
      IF(.NOT.ASSOCIATED(basis)) THEN
        CALL flagerror("Basis is not associated",err,error,*999)
      ENDIF

      SELECT CASE( basis%TYPE )
      CASE( basis_lagrange_hermite_tp_type )
        CALL field_io_calculate_tp_scale_and_node_counts(basis, num_scl, num_node, err, error, *999 )
      CASE( basis_simplex_type )
        CALL fieldio_calculatesimplexscaleandnodecounts(basis, num_scl, num_node, err, error, *999 )
      CASE DEFAULT
        CALL flagerror("Basis type "//trim(number_to_vstring(basis%TYPE,"*",err,error))//" is invalid or not implemented",&
            &err,error,*999)
      END SELECT

      group_scale_factors(scaleindex)=num_scl !numer of scale factors in scale factor set
      group_node(scaleindex)=num_node !numer of nodes in scale factor set
    ENDDO !scaleIndex

    !write out the scale factor set information
    err = fieldexport_scalingfactorcount( sessionhandle, num_of_scaling_factor_sets )
    IF(err/=0) THEN
      CALL flagerror( "File write error during field export", err, error,*999 )
    ENDIF

    CALL reallocate(interpolation_xi, basis%NUMBER_OF_XI, &
        & "Could not allocate temporary variable buffer in IO", err, error, *999)
    CALL reallocate(element_derivatives, sum(group_scale_factors(:)), &
        & "Could not allocate temporary variable buffer in IO", err, error, *999)


    DO scaleindex = 1, num_of_scaling_factor_sets
      basis => listscalebases( scaleindex )%PTR
      SELECT CASE( basis%TYPE )
      CASE( basis_lagrange_hermite_tp_type, basis_simplex_type )
        !!TEMP
        !ERR = FieldExport_ScaleFactors( sessionHandle, basis%NUMBER_OF_XI, C_LOC(basis%INTERPOLATION_XI) );
        !!Copy interpolation xi to a temporary array that has the target attribute. gcc bug 38813 prevents using C_LOC with
        !!the array directly. nb using a fixed length array here which is dangerous but should suffice for now.
        interpolation_xi(1:basis%NUMBER_OF_XI)=basis%INTERPOLATION_XI(1:basis%NUMBER_OF_XI)
        err = fieldexport_scalefactors( sessionhandle, basis%NUMBER_OF_XI, c_loc(interpolation_xi), &
          & basis%NUMBER_OF_ELEMENT_PARAMETERS );
        IF( err /= 0 ) THEN
          CALL flagerror( "can not get basis type of lagrange_hermite label" ,err, error, *999 )
        ENDIF
      CASE DEFAULT
        CALL flagerror( "Basis type "//trim(number_to_vstring(basis%TYPE, "*" , err, error ))//" is not implemented",&
            &err,error, *999)
      END SELECT
    ENDDO !scaleIndex

    err = fieldexport_nodecount( sessionhandle, max_num_nodes )
    IF(err/=0) THEN
      CALL flagerror( "File write error during field export", err, error,*999 )
    ENDIF

    err = fieldexport_fieldcount( sessionhandle, num_of_variables )
    IF(err/=0) THEN
      CALL flagerror( "File write error during field export", err, error,*999 )
    ENDIF

    !write out the nodal header
    var_idx=0
    NULLIFY(variable_ptr)
    DO comp_idx=1,elementalinfoset%NUMBER_OF_COMPONENTS
      component => elementalinfoset%COMPONENTS(comp_idx)%PTR

      !grouping field variables and components together
      IF(.NOT.ASSOCIATED(variable_ptr,TARGET=component%FIELD_VARIABLE)) THEN !different variables
        var_idx=var_idx+1
        variable_ptr=>component%FIELD_VARIABLE
        !write out the field information

        fvar_name = char(variable_ptr%variable_label)
        length=len_trim(fvar_name)
        DO i=1,length
          cvar_name(i)=fvar_name(i:i)
        ENDDO !i
        cvar_name(length+1)=c_null_char

        IF( variable_ptr%FIELD%TYPE == field_geometric_type .AND. &
            & variable_ptr%VARIABLE_TYPE == field_u_variable_type ) THEN
          NULLIFY(coordinate_system)
          CALL field_coordinate_system_get(variable_ptr%FIELD,coordinate_system,err,error,*999)
          err = fieldexport_coordinatevariable( sessionhandle, cvar_name, var_idx, coordinate_system%TYPE, &
              & group_variables(var_idx) )
        ELSE
          err = fieldexport_variable( sessionhandle, cvar_name, var_idx, variable_ptr%FIELD%TYPE, variable_ptr%VARIABLE_TYPE, &
              & group_variables(var_idx) )
        ENDIF

        IF( err /= 0 ) THEN
          CALL flagerror( "File write error during field export", err, error,*999 )
        ENDIF
      ENDIF

      componentdomain=>component%DOMAIN
      domain_elements=>componentdomain%TOPOLOGY%ELEMENTS
      basis=>domain_elements%ELEMENTS(group_local_number(comp_idx))%BASIS

      SELECT CASE( basis%TYPE )
      CASE( basis_lagrange_hermite_tp_type )
        CALL field_io_calculate_tp_scale_and_node_counts(basis, num_scl, num_node, err, error, *999 )
      CASE( basis_simplex_type )
        CALL fieldio_calculatesimplexscaleandnodecounts(basis, num_scl, num_node, err, error, *999 )
      CASE DEFAULT
        CALL flagerror("Basis type "//trim(number_to_vstring(basis%TYPE,"*",err,error))//" is invalid or not implemented",&
            &err,error,*999)
      END SELECT
      CALL reallocate(number_of_derivatives, num_node, &
          & "Could not allocate temporary variable buffer in IO", err, error, *999)
      CALL reallocate(scale_indexes, num_scl, &
          & "Could not allocate temporary variable buffer in IO", err, error, *999)
      CALL reallocate(node_indexes, num_node, &
          & "Could not allocate temporary variable buffer in IO", err, error, *999)
      
      SELECT CASE(component%INTERPOLATION_TYPE)
      CASE(field_constant_interpolation)
        interptype = 1
      CASE(field_element_based_interpolation)
        interptype = 2
      CASE(field_node_based_interpolation)
        interptype = 3
      CASE(field_grid_point_based_interpolation)
        interptype = 4
      CASE(field_gauss_point_based_interpolation)
        interptype = 5
      CASE(field_data_point_based_interpolation)
        interptype = 6
      CASE DEFAULT
        interptype = 0
      END SELECT

      IF(component%INTERPOLATION_TYPE==field_gauss_point_based_interpolation) THEN
        !TEMP HACK. Fake gauss point export as regular grid. Use interpolation xi to pass in number of Gauss.
        interpolation_xi(1:basis%NUMBER_OF_XI)=basis%QUADRATURE%NUMBER_OF_GAUSS_XI(1:basis%NUMBER_OF_XI)
      ELSE
        !Copy interpolation xi to a temporary array that has the target attribute. gcc bug 38813 prevents using C_LOC with
        !the array directly. nb using a fixed length array here which is dangerous but should suffice for now.
        interpolation_xi(1:basis%NUMBER_OF_XI)=basis%INTERPOLATION_XI(1:basis%NUMBER_OF_XI)
      ENDIF
      
      IF( variable_ptr%FIELD%TYPE == field_geometric_type .AND. &
          & variable_ptr%VARIABLE_TYPE == field_u_variable_type ) THEN
        !!TEMP
        !ERR = FieldExport_CoordinateComponent( sessionHandle, variable_ptr%FIELD%REGION%COORDINATE_SYSTEM, &
        !  & component%COMPONENT_NUMBER, basis%NUMBER_OF_XI, C_LOC( basis%INTERPOLATION_XI ) )
        NULLIFY(coordinate_system)
        CALL field_coordinate_system_get(variable_ptr%FIELD,coordinate_system,err,error,*999)
        err = fieldexport_coordinatecomponent( sessionhandle, coordinate_system%TYPE, &
            & component%COMPONENT_NUMBER,interptype,basis%NUMBER_OF_XI, c_loc( interpolation_xi ))
      ELSE
        !!TEMP
        !ERR = FieldExport_Component( sessionHandle, &
        !  & component%COMPONENT_NUMBER, basis%NUMBER_OF_XI, C_LOC( basis%INTERPOLATION_XI ) )
        err = fieldexport_component( sessionhandle, &
            & component%COMPONENT_NUMBER,interptype,basis%NUMBER_OF_XI, c_loc( interpolation_xi ) )
      ENDIF
      IF(err/=0) THEN
        CALL flagerror( "File write error during field export", err, error,*999 )
      ENDIF

      IF( interptype /= 3 .AND. interptype /= 6) THEN
        err = fieldexport_elementgridsize( sessionhandle, interptype, basis%NUMBER_OF_XI, c_loc( interpolation_xi ) )
      ELSE
        ! IF(.NOT.BASIS%DEGENERATE) THEN
        IF(list_comp_scale(comp_idx)==1) THEN
          scaleindex=0
        ELSE
          scaleindex= sum(group_scale_factors(1:list_comp_scale(comp_idx)))-1
        ENDIF

        ! Fortran numbering instead of c numbering
        scaleindex1 = scaleindex + 1

        !!TEMP
        ! ERR = FieldExport_NodeScaleIndexes( sessionHandle, BASIS%NUMBER_OF_NODES, C_LOC( BASIS%NUMBER_OF_DERIVATIVES ), &
        ! & C_LOC( DOMAIN_ELEMENTS%ELEMENTS(GROUP_LOCAL_NUMBER(comp_idx))%ELEMENT_DERIVATIVES ), scaleIndex )
        !!Copy element derivatives etc. to a temporary array that has the target attribute. gcc bug 38813 prevents using C_LOC with
        !!the array directly. nb using a fixed length array here which is dangerous but should suffice for now.
        !!In order to correctly index the supplied array, the API needs to know in advance the dimensions of the array.
        !!To avoid having to pass in an extra 'size' parameter, we unroll the 2d derivative index array into a vector.
        derivativeindex = 1
        nn=0
        node_number=0

        ! For elements with collapsed nodes, the node indexes need to be changed
        IF (basis%NUMBER_OF_COLLAPSED_XI>0) THEN
    !!$         IF((.NOT.BASIS%COLLAPSED_XI(1).OR.BASIS%COLLAPSED_XI(2).OR.BASIS%COLLAPSED_XI(3)==BASIS_NOT_COLLAPSED)) THEN
    !!$           CALL FlagError("Pyramide elements output is not implemented",&
    !!$              &ERR,ERROR,*999)

          node_number_counter=0
          number_of_unique_nodes = 0
          IF(basis%INTERPOLATION_XI(1)>3) THEN
            nodesx=2
          ELSE
            nodesx=basis%INTERPOLATION_XI(1)+1
          ENDIF
          IF(basis%INTERPOLATION_XI(2)>3) THEN
            nodesy=2
          ELSE
            nodesy=basis%INTERPOLATION_XI(2)+1
          ENDIF
          IF(basis%INTERPOLATION_XI(3)>3) THEN
            nodesz=2
          ELSE
            nodesz=basis%INTERPOLATION_XI(3)+1
          ENDIF

          !The following if-sentences goes through all possible wedge formed elements and renumber the nodes in order to
          !attach the node_index and the number of derivatives to the numbering corresponding to not collapsed elements
          IF(basis%COLLAPSED_XI(1)==basis_xi_collapsed) THEN
            IF(basis%COLLAPSED_XI(2)==basis_collapsed_at_xi0) THEN
              DO nz=1,nodesz
                DO ny=1,nodesy
                  DO nx=1,nodesx
                    nn=nn+1
                    IF (ny==1) THEN
                      node_number_collapsed=(nz-1)*nodesx*(nodesy-1)+nz
                      node_number=node_number_collapsed
                    ELSE
                      IF (node_number_counter<node_number_collapsed) THEN
                        node_number_counter=node_number_collapsed+1
                      ELSE
                        node_number_counter=node_number_counter+1
                      ENDIF
                      node_number=node_number_counter
                    ENDIF
                    node_indexes(nn)=node_number
                    number_of_derivatives(nn)=basis%NUMBER_OF_DERIVATIVES(node_number)
                    IF(node_number>number_of_unique_nodes) THEN
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex)=scaleindex1 
                        derivativeindex=derivativeindex+1
                        scaleindex1=scaleindex1+1
                      ENDDO
                      number_of_unique_nodes=number_of_unique_nodes+1
                    ELSE
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex) = scale_indexes(sum(number_of_derivatives(1:node_number+ &
                          & nn-number_of_unique_nodes-2))+mm)
                        derivativeindex = derivativeindex + 1
                      ENDDO
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
            ELSE IF(basis%COLLAPSED_XI(2)==basis_collapsed_at_xi1) THEN
              DO nz=1,nodesz
                DO ny=1,nodesy
                  DO nx=1,nodesx
                    nn=nn+1
                    IF (ny==nodesy) THEN
                      node_number_collapsed=(nz-1)*nodesx*(nodesy-1)+nodesx*(nodesy-1)+nz
                      node_number=node_number_collapsed
                    ELSE
                      IF (node_number_counter<node_number_collapsed) THEN
                        node_number_counter=node_number_collapsed+1
                      ELSE
                        node_number_counter=node_number_counter+1
                      ENDIF
                      node_number=node_number_counter
                    ENDIF
                    node_indexes(nn)=node_number
                    number_of_derivatives(nn)=basis%NUMBER_OF_DERIVATIVES(node_number)
                    IF(node_number>number_of_unique_nodes) THEN
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex)=scaleindex1 
                        derivativeindex=derivativeindex+1
                        scaleindex1=scaleindex1+1
                      ENDDO
                      number_of_unique_nodes=number_of_unique_nodes+1
                    ELSE
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex) = scale_indexes(sum(number_of_derivatives(1:node_number+ &
                          & nn-number_of_unique_nodes-2))+mm)
                        derivativeindex = derivativeindex + 1
                      ENDDO
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
            ELSE IF(basis%COLLAPSED_XI(3)==basis_collapsed_at_xi0) THEN
              DO nz=1,nodesz
                DO ny=1,nodesy
                  DO nx=1,nodesx
                    nn=nn+1
                    IF (nz==1) THEN
                      node_number_collapsed=ny
                      node_number=node_number_collapsed
                    ELSE
                      IF (node_number_counter<node_number_collapsed) THEN
                        node_number_counter=node_number_collapsed+1
                      ELSE
                        node_number_counter=node_number_counter+1
                      ENDIF
                      node_number=node_number_counter
                    ENDIF
                    node_indexes(nn)=node_number
                    number_of_derivatives(nn)=basis%NUMBER_OF_DERIVATIVES(node_number)
                    IF(node_number>number_of_unique_nodes) THEN
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex)=scaleindex1 
                        derivativeindex=derivativeindex+1
                        scaleindex1=scaleindex1+1
                      ENDDO
                      number_of_unique_nodes=number_of_unique_nodes+1
                    ELSE
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex) = scale_indexes(sum(number_of_derivatives(1:node_number+ &
                          & nn-number_of_unique_nodes-2))+mm)
                        derivativeindex = derivativeindex + 1
                      ENDDO
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
            ELSE IF(basis%COLLAPSED_XI(3)==basis_collapsed_at_xi1) THEN
              DO nz=1,nodesz
                DO ny=1,nodesy
                  DO nx=1,nodesx
                    nn=nn+1
                    IF (nz==nodesz) THEN
                      node_number_collapsed=(nodesz-1)*(nodesy*nodesx)+ny
                      node_number=node_number_collapsed
                    ELSE
                      IF (node_number_counter<node_number_collapsed) THEN
                        node_number_counter=node_number_collapsed+1
                      ELSE
                        node_number_counter=node_number_counter+1
                      ENDIF
                      node_number=node_number_counter
                    ENDIF
                    node_indexes(nn)=node_number
                    number_of_derivatives(nn)=basis%NUMBER_OF_DERIVATIVES(node_number)
                    IF(node_number>number_of_unique_nodes) THEN
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex)=scaleindex1 
                        derivativeindex=derivativeindex+1
                        scaleindex1=scaleindex1+1
                      ENDDO
                      number_of_unique_nodes=number_of_unique_nodes+1
                    ELSE
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex) = scale_indexes(sum(number_of_derivatives(1:node_number+ &
                          & nn-number_of_unique_nodes-2))+mm)
                        derivativeindex = derivativeindex + 1
                      ENDDO
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
            ENDIF
          ELSE IF(basis%COLLAPSED_XI(2)==basis_xi_collapsed) THEN
            IF(basis%COLLAPSED_XI(1)==basis_collapsed_at_xi0) THEN
              DO nz=1,nodesz
                DO ny=1,nodesy
                  DO nx=1,nodesx
                    nn=nn+1
                    IF (ny==1) THEN
                      node_number_collapsed=(nz-1)*nodesx*(nodesy-1)+nz
                      node_number=node_number_collapsed
                    ELSE
                      IF (node_number_counter<node_number_collapsed) THEN
                        node_number_counter=node_number_collapsed+1
                      ELSE
                        node_number_counter=node_number_counter+1
                      ENDIF
                      node_number=node_number_counter
                    ENDIF
                    node_indexes(nn)=node_number
                    number_of_derivatives(nn)=basis%NUMBER_OF_DERIVATIVES(node_number)
                    IF(node_number>number_of_unique_nodes) THEN
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex)=scaleindex1 
                        derivativeindex=derivativeindex+1
                        scaleindex1=scaleindex1+1
                      ENDDO
                      number_of_unique_nodes=number_of_unique_nodes+1
                    ELSE
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex) = scale_indexes(sum(number_of_derivatives(1:node_number+ &
                          & nn-number_of_unique_nodes-2))+mm)
                        derivativeindex = derivativeindex + 1
                      ENDDO
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
            ELSE IF(basis%COLLAPSED_XI(1)==basis_collapsed_at_xi1) THEN
              DO nz=1,nodesz
                DO ny=1,nodesy
                  DO nx=1,nodesx
                    nn=nn+1
                    IF (ny==nodesy) THEN
                      node_number=(nz-1)*nodesx*(nodesy-1)+nz+nodesx-1
                      node_number=node_number_collapsed
                    ELSE
                      IF (node_number_counter<node_number_collapsed) THEN
                        node_number_counter=node_number_collapsed+1
                      ELSE
                        node_number_counter=node_number_counter+1
                      ENDIF
                      node_number=node_number_counter
                    ENDIF
                    node_indexes(nn)=node_number
                    number_of_derivatives(nn)=basis%NUMBER_OF_DERIVATIVES(node_number)
                    IF(node_number>number_of_unique_nodes) THEN
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex)=scaleindex1 
                        derivativeindex=derivativeindex+1
                        scaleindex1=scaleindex1+1
                      ENDDO
                      number_of_unique_nodes=number_of_unique_nodes+1
                    ELSE
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex) = scale_indexes(sum(number_of_derivatives(1:node_number+ &
                          & nn-number_of_unique_nodes-2))+mm)
                        derivativeindex = derivativeindex + 1
                      ENDDO
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
            ELSE IF(basis%COLLAPSED_XI(3)==basis_collapsed_at_xi0) THEN
              DO nz=1,nodesz
                DO ny=1,nodesy
                  DO nx=1,nodesx
                    nn=nn+1
                    IF (nz==1) THEN
                      node_number_collapsed=nx
                      node_number=node_number_collapsed
                    ELSE
                      IF (node_number_counter<node_number_collapsed) THEN
                        node_number_counter=node_number_collapsed+1
                      ELSE
                        node_number_counter=node_number_counter+1
                      ENDIF
                      node_number=node_number_counter
                    ENDIF
                    node_indexes(nn)=node_number
                    number_of_derivatives(nn)=basis%NUMBER_OF_DERIVATIVES(node_number)
                    IF(node_number>number_of_unique_nodes) THEN
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex)=scaleindex1 
                        derivativeindex=derivativeindex+1
                        scaleindex1=scaleindex1+1
                      ENDDO
                      number_of_unique_nodes=number_of_unique_nodes+1
                    ELSE
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex) = scale_indexes(sum(number_of_derivatives(1:node_number+ &
                          & nn-number_of_unique_nodes-2))+mm)
                        derivativeindex = derivativeindex + 1
                      ENDDO
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
            ELSE IF(basis%COLLAPSED_XI(3)==basis_collapsed_at_xi1) THEN
              DO nz=1,nodesz
                DO ny=1,nodesy
                  DO nx=1,nodesx
                    nn=nn+1
                    IF (nz==nodesz) THEN
                      node_number_collapsed=(nodesz-1)*nodesy*nodesx+nx
                      node_number=node_number_collapsed
                    ELSE
                      IF (node_number_counter<node_number_collapsed) THEN
                        node_number_counter=node_number_collapsed+1
                      ELSE
                        node_number_counter=node_number_counter+1
                      ENDIF
                      node_number=node_number_counter
                    ENDIF
                    node_indexes(nn)=node_number
                    number_of_derivatives(nn)=basis%NUMBER_OF_DERIVATIVES(node_number)
                    IF(node_number>number_of_unique_nodes) THEN
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex)=scaleindex1 
                        derivativeindex=derivativeindex+1
                        scaleindex1=scaleindex1+1
                      ENDDO
                      number_of_unique_nodes=number_of_unique_nodes+1
                    ELSE
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex) = scale_indexes(sum(number_of_derivatives(1:node_number+ &
                          & nn-number_of_unique_nodes-2))+mm)
                        derivativeindex = derivativeindex + 1
                      ENDDO
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
            ENDIF
          ELSE IF(basis%COLLAPSED_XI(3)==basis_xi_collapsed) THEN
            IF(basis%COLLAPSED_XI(1)==basis_collapsed_at_xi0) THEN
              DO nz=1,nodesz
                DO ny=1,nodesy
                  DO nx=1,nodesx
                    nn=nn+1
                    IF (nx==1) THEN
                      node_number_collapsed=nodesx*(ny-1)+ny
                      node_number=node_number_collapsed
                    ELSE
                      IF (node_number_counter<node_number_collapsed) THEN
                        node_number_counter=node_number_collapsed+1
                      ELSE
                        node_number_counter=node_number_counter+1
                      ENDIF
                      node_number=node_number_counter
                    ENDIF
                    node_indexes(nn)=node_number
                    number_of_derivatives(nn)=basis%NUMBER_OF_DERIVATIVES(node_number)
                    IF(node_number>number_of_unique_nodes) THEN
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex)=scaleindex1 
                        derivativeindex=derivativeindex+1
                        scaleindex1=scaleindex1+1
                      ENDDO
                      number_of_unique_nodes=number_of_unique_nodes+1
                    ELSE
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex) = scale_indexes(sum(number_of_derivatives(1:node_number+ &
                          & nn-number_of_unique_nodes-2))+mm)
                        derivativeindex = derivativeindex + 1
                      ENDDO
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
            ELSE IF(basis%COLLAPSED_XI(1)==basis_collapsed_at_xi1) THEN
              DO nz=1,nodesz
                DO ny=1,nodesy
                  DO nx=1,nodesx
                    nn=nn+1
                    IF (nx==nodesx) THEN
                      node_number_collapsed=nodesx*(ny-1)+ny+nodesx-1
                      node_number=node_number_collapsed
                    ELSE
                      IF (node_number_counter<node_number_collapsed) THEN
                        node_number_counter=node_number_collapsed+1
                      ELSE
                        node_number_counter=node_number_counter+1
                      ENDIF
                      node_number=node_number_counter
                    ENDIF
                    node_indexes(nn)=node_number
                    number_of_derivatives(nn)=basis%NUMBER_OF_DERIVATIVES(node_number)
                    IF(node_number>number_of_unique_nodes) THEN
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex)=scaleindex1 
                        derivativeindex=derivativeindex+1
                        scaleindex1=scaleindex1+1
                      ENDDO
                      number_of_unique_nodes=number_of_unique_nodes+1
                    ELSE
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex) = scale_indexes(sum(number_of_derivatives(1:node_number+ &
                          & nn-number_of_unique_nodes-2))+mm)
                        derivativeindex = derivativeindex + 1
                      ENDDO
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
            ELSE IF(basis%COLLAPSED_XI(2)==basis_collapsed_at_xi0) THEN
              DO nz=1,nodesz
                DO ny=1,nodesy
                  DO nx=1,nodesx
                    nn=nn+1
                    IF (ny==1) THEN
                      node_number_collapsed=nx
                      node_number=node_number_collapsed
                    ELSE
                      IF (node_number_counter<node_number_collapsed) THEN
                        node_number_counter=node_number_collapsed+1
                      ELSE
                        node_number_counter=node_number_counter+1
                      ENDIF
                      node_number=node_number_counter
                    ENDIF
                    node_indexes(nn)=node_number
                    number_of_derivatives(nn)=basis%NUMBER_OF_DERIVATIVES(node_number)
                    IF(node_number>number_of_unique_nodes) THEN
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex)=scaleindex1 
                        derivativeindex=derivativeindex+1
                        scaleindex1=scaleindex1+1
                      ENDDO
                      number_of_unique_nodes=number_of_unique_nodes+1
                    ELSE
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex) = scale_indexes(sum(number_of_derivatives(1:node_number+ &
                          & nn-number_of_unique_nodes-2))+mm)
                        derivativeindex = derivativeindex + 1
                      ENDDO
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
            ELSE IF(basis%COLLAPSED_XI(2)==basis_collapsed_at_xi1) THEN
              DO nz=1,nodesz
                DO ny=1,nodesy
                  DO nx=1,nodesx
                    nn=nn+1
                    IF (ny==nodesy) THEN
                      node_number_collapsed=nodesx*(nodesy-1)+nx
                      node_number=node_number_collapsed
                    ELSE
                      IF (node_number_counter<node_number_collapsed) THEN
                        node_number_counter=node_number_collapsed+1
                      ELSE
                        node_number_counter=node_number_counter+1
                      ENDIF
                      node_number=node_number_counter
                    ENDIF
                    node_indexes(nn)=node_number
                    number_of_derivatives(nn)=basis%NUMBER_OF_DERIVATIVES(node_number)
                    IF(node_number>number_of_unique_nodes) THEN
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex)=scaleindex1 
                        derivativeindex=derivativeindex+1
                        scaleindex1=scaleindex1+1
                      ENDDO
                      number_of_unique_nodes=number_of_unique_nodes+1
                    ELSE
                      DO mm=1,number_of_derivatives(nn)
                        element_derivatives(derivativeindex)=field_io_element_derivative_index( &
                          & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,node_number,err,error)
                        scale_indexes(derivativeindex) = scale_indexes(sum(number_of_derivatives(1:node_number+ &
                          & nn-number_of_unique_nodes-2))+mm)
                        derivativeindex = derivativeindex + 1
                      ENDDO
                    ENDIF
                  ENDDO
                ENDDO
              ENDDO
            ENDIF
          ENDIF
          !NODE_NUMBER_COUNTER has counted up to the highest number of nodes in the element
          number_of_element_nodes=nn
          !  ENDIF
        ELSE !Not collapsed nodes
          DO nn=1,basis%NUMBER_OF_NODES
            number_of_derivatives(nn) = basis%NUMBER_OF_DERIVATIVES(nn)
            DO mm=1,number_of_derivatives(nn)
              element_derivatives(derivativeindex) = field_io_element_derivative_index( &
                & domain_elements%ELEMENTS(group_local_number(comp_idx)),mm,nn,err,error)
              scale_indexes(derivativeindex) = scaleindex1
              derivativeindex = derivativeindex + 1
              scaleindex1 = scaleindex1 + 1
            ENDDO !mm
          ENDDO !nn

          !Find the local-node index in the element's total node list.
          !TODO This assumes nested subsets of nodes, and will therefore break on, e.g., mixed quad and cubic interpolation
          DO nn = 1, basis%NUMBER_OF_NODES
            DO mm = 1, max_node_element%BASIS%NUMBER_OF_NODES
              node_local_number = domain_elements%ELEMENTS( local_number )%ELEMENT_NODES( nn )
              node_user_number=domain_elements%DOMAIN%TOPOLOGY%NODES%NODES(node_local_number)%USER_NUMBER
              max_element_local_number = max_node_element%ELEMENT_NODES( mm )
              max_element_user_number = max_element_domain_nodes%NODES(max_element_local_number)%USER_NUMBER
              IF( node_user_number == max_element_user_number ) THEN
                node_indexes( nn ) = mm
                EXIT
              ENDIF
            ENDDO !mm
          ENDDO !nn
          number_of_element_nodes= basis%NUMBER_OF_NODES
        ENDIF


        IF( variable_ptr%FIELD%SCALINGS%SCALING_TYPE == field_no_scaling ) THEN
          scale_indexes(:) = -1
        ENDIF
        err = fieldexport_nodescaleindexes( sessionhandle, number_of_element_nodes, c_loc( number_of_derivatives ), &
            & c_loc( element_derivatives ), c_loc( node_indexes ), c_loc( scale_indexes ) )
        ! ELSE
        !  CALL FlagError("exporting degenerated nodes has not been implemented",ERR,ERROR,*999)
        ! ENDIF
      ENDIF

      IF(err/=0) THEN
        CALL flagerror( "File write error during field export", err, error,*999 )
      ENDIF

    ENDDO !comp_idx

    !release temporary memory
    CALL checked_deallocate( listscalefields )
    CALL checked_deallocate( listscalebases )
    CALL checked_deallocate( group_local_number )
    CALL checked_deallocate( group_scale_factors )
    CALL checked_deallocate( group_node )
    CALL checked_deallocate( group_variables )
    CALL checked_deallocate( interpolation_xi )
    CALL checked_deallocate( element_derivatives )
    CALL checked_deallocate( number_of_derivatives )
    CALL checked_deallocate( node_indexes )
    CALL checked_deallocate( scale_indexes )


    exits("FieldIO_ExportElementalGroupHeaderFortran")
    RETURN
999 errorsexits("FieldIO_ExportElementalGroupHeaderFortran",err,error)
    RETURN 1
  END SUBROUTINE fieldio_exportelementalgroupheaderfortran

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

  SUBROUTINE field_io_export_element_scale_factors( sessionHandle, components, componentScales, globalNumber, &
    & mycomputationalnodenumber, err, error, * )
    !Argument variables
    INTEGER(INTG) :: sessionHandle
    TYPE(field_io_component_info_set), INTENT(INOUT) :: components
    INTEGER(INTG) :: componentScales(:)
    INTEGER(INTG) :: globalNumber
    INTEGER(INTG) :: myComputationalNodeNumber
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    !Local variables
    INTEGER(INTG) :: scaleIndex, componentIndex, localNumber, scaleFactorCount, nodeIndex
    INTEGER(INTG) :: nodeNumber, derivativeIndex, nv, nk, ny2, firstScaleSet
    TYPE(field_variable_component_type), POINTER :: component
    TYPE(domain_elements_type), POINTER :: domainElements
    TYPE(domain_nodes_type), POINTER :: domainNodes
    TYPE(domain_mapping_type), POINTER :: domainElementMapping
    TYPE(basis_type), POINTER :: basis
    REAL(C_DOUBLE), ALLOCATABLE, TARGET :: scaleBuffer(:)
    REAL(DP), POINTER :: SCALE_FACTORS(:)

    NULLIFY(scale_factors)

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

    scaleindex = 1
    firstscaleset = 1
    DO componentindex = 1, components%NUMBER_OF_COMPONENTS
      component => components%COMPONENTS( componentindex )%PTR
      !finding the local numbering through the global to local mapping

      domainelementmapping=>component%DOMAIN%MAPPINGS%ELEMENTS
      !get the domain index for this variable component according to my own computional node number

      localnumber = findmylocaldomainnumber( domainelementmapping%GLOBAL_TO_LOCAL_MAP( globalnumber ), &
        & mycomputationalnodenumber )
      !use local domain information find the out the maximum number of derivatives
      domainelements => component%DOMAIN%TOPOLOGY%ELEMENTS
      domainnodes => component%DOMAIN%TOPOLOGY%NODES

      !write out the components' values of this node in this domain
      !DO scaleIndex=1, NUM_OF_SCALING_FACTOR_SETS
      IF( componentscales( componentindex ) == scaleindex ) THEN
        scaleindex = scaleindex + 1

        scalefactorcount = 0
        basis => domainelements%ELEMENTS( localnumber )%BASIS

        CALL reallocate( scalebuffer, sum( basis%NUMBER_OF_DERIVATIVES(1:basis%NUMBER_OF_NODES ) ), &
          & "Could not allocate scale buffer in IO", err, error, *999 )

        IF( component%FIELD_VARIABLE%FIELD%SCALINGS%SCALING_TYPE /= field_no_scaling ) THEN
          CALL distributed_vector_data_get(component%FIELD_VARIABLE%FIELD%SCALINGS%SCALINGS(component% &
            & scaling_index)%SCALE_FACTORS,scale_factors,err,error,*999)
        ENDIF

        !IF( .NOT.basis%DEGENERATE ) THEN
          DO nodeindex = 1, basis%NUMBER_OF_NODES
            nodenumber = domainelements%ELEMENTS( localnumber )%ELEMENT_NODES( nodeindex )
            DO derivativeindex = 1, basis%NUMBER_OF_DERIVATIVES( nodeindex )
              nk = domainelements%ELEMENTS( localnumber )%ELEMENT_DERIVATIVES(derivativeindex, nodeindex )
              nv = domainelements%ELEMENTS( localnumber )%elementVersions(derivativeindex, nodeindex )
              ny2 = domainnodes%NODES( nodenumber )%DERIVATIVES(nk)%DOF_INDEX(nv)
              scalefactorcount = scalefactorcount + 1
              IF( component%FIELD_VARIABLE%FIELD%SCALINGS%SCALING_TYPE /= field_no_scaling ) THEN
                scalebuffer( scalefactorcount ) = scale_factors(ny2)
              ELSE
                scalebuffer( scalefactorcount ) = 1
              ENDIF
            ENDDO !derivativeIndex
          ENDDO !nodeIndex
        !ELSE
          !This is just a hack, forcing to write out the correct number of scale factors equal to one!!!!
        !  NodesX=BASIS%INTERPOLATION_XI(1)+1
        !  NodesY=BASIS%INTERPOLATION_XI(2)+1
        !  NodesZ=BASIS%INTERPOLATION_XI(3)+1
        !  CALL REALLOCATE( scaleBuffer, (NodesX*NodesY*NodesZ), &
        !       & "Could not allocate scale buffer in IO", ERR, ERROR, *999 )
        !  DO nz=1,NodesZ
        !    DO ny=1,NodesY
        !      DO nx=1,NodesX
        !        scaleFactorCount=scaleFactorCOUNT+1
        !        scaleBuffer( scaleFactorCount ) = 1
        !      ENDDO
        !    ENDDO
        !  ENDDO
        !ENDIF

        NULLIFY( scale_factors )

        err = fieldexport_elementnodescales( sessionhandle, firstscaleset, scalefactorcount, c_loc( scalebuffer ) )

        firstscaleset = 0

        IF( err /= 0 ) THEN
          CALL flagerror( "Cannot write node scales to file", err, error,*999 )
        ENDIF

      ENDIF ! componentScales(componentIndex) == scaleIndex
    ENDDO ! componentIndex

    CALL checked_deallocate( scalebuffer )

    exits("FIELD_IO_EXPORT_ELEMENT_SCALE_FACTORS")
    RETURN
999 errorsexits("FIELD_IO_EXPORT_ELEMENT_SCALE_FACTORS",err,error)
    RETURN 1
  END SUBROUTINE field_io_export_element_scale_factors

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

  SUBROUTINE field_io_export_elements_into_local_file(ELEMENTAL_INFO_SET, NAME, my_computational_node_number, &
      & err, error, *)
    !the reason that my_computational_node_number is used in the argument is for future extension
    !Argument variables
    TYPE(field_io_info_set), INTENT(INOUT) :: ELEMENTAL_INFO_SET
    TYPE(varying_string), INTENT(IN) :: NAME
    INTEGER(INTG), INTENT(IN):: my_computational_node_number
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: sessionHandle
    TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
    TYPE(field_variable_component_type), POINTER :: component
    TYPE(mesh_element_type), POINTER :: element
    TYPE(varying_string) :: FILE_NAME !the prefix name of file.
    TYPE(basis_type), POINTER ::BASIS
    TYPE(domain_mapping_type), POINTER :: DOMAIN_MAPPING_ELEMENTS !The domain mapping to calculate elemental mappings
    TYPE(domain_elements_type), POINTER :: DOMAIN_ELEMENTS ! domain elements
    INTEGER(INTG) :: local_number, global_number, MAX_NODE_COMP_INDEX, NUM_DIM
    INTEGER(INTG), ALLOCATABLE :: LIST_COMP_SCALE(:), NODAL_NUMBER(:)!LIST_COMP(:) !Components which will be used for export scale factors
    INTEGER(C_INT), TARGET :: USER_ELEMENT_NODES(64)
    INTEGER(INTG) :: elem_idx, comp_idx, NUM_OF_SCALING_FACTOR_SETS, isFirstValueSet !dev_idx  elem_num
    REAL(DP), ALLOCATABLE :: SCALE_FACTORS(:)
    TYPE(field_io_component_info_set), POINTER :: components
    REAL(DP), POINTER :: GEOMETRIC_PARAMETERS(:)
    INTEGER(INTG), POINTER :: GEOMETRIC_PARAMETERS_INTG(:)
    REAL(DP), ALLOCATABLE :: GEOMETRIC_PARAMETERS_DP(:)

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

    !is not necessarily equal to numbering of computional node, so use method COMPUTATIONAL_NODE_NUMBER_GET
    !will be a secured way to get the number
    !my_computational_node_number=COMPUTATIONAL_NODE_NUMBER_GET(ERR,ERROR)
    !IF(ERR/=0) GOTO 999
    file_name=name//".part"//trim(number_to_vstring(my_computational_node_number,"*",err,error))//".exelem"
    num_of_scaling_factor_sets=0

    IF(.NOT.ALLOCATED(elemental_info_set%COMPONENT_INFO_SET)) THEN
      CALL flagerror("the elemental information set in input is invalid",err,error,*999)
    ENDIF

    IF(.NOT.ALLOCATED(elemental_info_set%LIST_OF_GLOBAL_NUMBER)) THEN
      CALL flagerror("the elemental information set is not associated with any numbering list",err,error,*999)
    ENDIF

    IF(elemental_info_set%NUMBER_OF_ENTRIES==0) THEN
      CALL flagerror("the elemental information set does not contain any nodes",err,error,*999)
    ENDIF

    IF(elemental_info_set%COMPONENT_INFO_SET(1)%PTR%SAME_HEADER) THEN
      CALL flagerror("the first header flag of elemental information set should be false",err,error,*999)
    ENDIF

    !NULLIFY(SCALE_FACTORS)
    !NULLIFY(LIST_COMP_SCALE)
    !NULLIFY(tmp_components)

    NULLIFY(coordinate_system)
    CALL field_coordinate_system_get(elemental_info_set%COMPONENT_INFO_SET(1)%PTR%COMPONENTS(1)%PTR% &
        & field_variable%FIELD,coordinate_system,err,error,*999)
    num_dim=coordinate_system%NUMBER_OF_DIMENSIONS

    err = fieldexport_opensession( export_type_file, char(file_name)//c_null_char, sessionhandle )
    IF(err/=0) THEN
      CALL flagerror( "Cannot open file export session", err, error,*999 )
    ENDIF

    IF(ASSOCIATED(elemental_info_set%FIELDS%REGION)) THEN
      err = fieldexport_group( sessionhandle, char(elemental_info_set%FIELDS%REGION%LABEL)//c_null_char )
    ELSE
      IF(ASSOCIATED(elemental_info_set%FIELDS%INTERFACE)) THEN
        err = fieldexport_group( sessionhandle, char(elemental_info_set%FIELDS%INTERFACE%LABEL)//c_null_char )
      ELSE
        CALL flagerror("Fields region or interface is not associated.",err,error,*999)
      ENDIF
    ENDIF
    IF(err/=0) THEN
      CALL flagerror( "Cannot write group name to elements file", err, error,*999 )
    ENDIF

    components => elemental_info_set%COMPONENT_INFO_SET(1)%PTR
    component => components%COMPONENTS(1)%PTR
    domain_mapping_elements=>component%DOMAIN%MAPPINGS%ELEMENTS
    domain_elements=>component%DOMAIN%TOPOLOGY%ELEMENTS
    basis => domain_elements%ELEMENTS( 1 )%BASIS

    err = fieldexport_meshdimensions( sessionhandle, basis%NUMBER_OF_XI, basis%TYPE  )
    IF(err/=0) THEN
      CALL flagerror( "Cannot write mesh dimensions to file", err, error,*999 )
    ENDIF

    DO elem_idx=1, elemental_info_set%NUMBER_OF_ENTRIES

      components => elemental_info_set%COMPONENT_INFO_SET(elem_idx)%PTR
      global_number = elemental_info_set%LIST_OF_GLOBAL_NUMBER(elem_idx)

      IF(.NOT.ALLOCATED(list_comp_scale)) THEN
        ALLOCATE(list_comp_scale(components%NUMBER_OF_COMPONENTS),stat=err)
        IF(err/=0) CALL flagerror("Could not allocate LIST_COMP_SCALE in exelem io",err,error,*999)
      ENDIF

      !check whether need to write out the nodal information header
      IF(.NOT.components%SAME_HEADER) THEN
        !write out the nodal header
        CALL fieldio_exportelementalgroupheaderfortran( global_number, max_node_comp_index, num_of_scaling_factor_sets, &
          & list_comp_scale, my_computational_node_number, components, sessionhandle, err, error, *999)
      ENDIF

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !write out elemental information
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !element info
      component => components%COMPONENTS(max_node_comp_index)%PTR
      element => component%DOMAIN%MESH%TOPOLOGY(component%MESH_COMPONENT_NUMBER)%PTR%ELEMENTS%ELEMENTS(global_number)

      err = fieldexport_elementindex( sessionhandle, num_dim, element%USER_NUMBER )
      IF(err/=0) THEN
        CALL flagerror( "Cannot write element index to file", err, error,*999 )
      ENDIF

      isfirstvalueset = 1
      DO comp_idx = 1, components%NUMBER_OF_COMPONENTS
        component => components%COMPONENTS(comp_idx)%PTR

        !finding the local numbering through the global to local mapping
        domain_mapping_elements=>component%DOMAIN%MAPPINGS%ELEMENTS
        domain_elements=>component%DOMAIN%TOPOLOGY%ELEMENTS
        !get the domain index for this variable component according to my own computional node number
        local_number = findmylocaldomainnumber( domain_mapping_elements%GLOBAL_TO_LOCAL_MAP( global_number ), &
          & my_computational_node_number )
        !use local domain information find the out the maximum number of derivatives
        basis => domain_elements%ELEMENTS( local_number )%BASIS

        IF( component%INTERPOLATION_TYPE == field_element_based_interpolation ) THEN
!         IF( .NOT.ASSOCIATED( component%DOMAIN%TOPOLOGY%ELEMENTS ) ) THEN
!           CYCLE
!         ENDIF
          IF(component%FIELD_VARIABLE%DATA_TYPE==field_dp_type) THEN
            NULLIFY(geometric_parameters)
            CALL field_parameter_set_data_get(component%FIELD_VARIABLE%FIELD,&
                & component%FIELD_VARIABLE%VARIABLE_TYPE,field_values_set_type,geometric_parameters,err,error,*999)
            err = fieldexport_elementgridvalues( sessionhandle, isfirstvalueset, 1, &
                & geometric_parameters(component%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(local_number)))
          ELSE IF(component%FIELD_VARIABLE%DATA_TYPE==field_intg_type) THEN
            NULLIFY(geometric_parameters_intg)
            CALL field_parameter_set_data_get(component%FIELD_VARIABLE%FIELD,&
                & component%FIELD_VARIABLE%VARIABLE_TYPE,field_values_set_type,geometric_parameters_intg,err,error,*999)
            ALLOCATE(geometric_parameters_dp(SIZE(geometric_parameters_intg)))
            IF(err/=0) CALL flagerror("Could not allocate geometric parameters dp", err, error,*999 )
            geometric_parameters_dp(1:SIZE(geometric_parameters_intg))= &
              & REAL(geometric_parameters_intg(1:size(geometric_parameters_intg)))
            err = fieldexport_elementgridvalues( sessionhandle, isfirstvalueset, 1, &
                & geometric_parameters_dp(component%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(local_number)))
            DEALLOCATE(geometric_parameters_dp)
          ELSE
            CALL flagerror( "Only INTG and REAL data types implemented.", err, error,*999 )
          ENDIF
          isfirstvalueset = 0
        ELSE IF( component%INTERPOLATION_TYPE == field_constant_interpolation ) THEN
          IF(component%FIELD_VARIABLE%DATA_TYPE==field_dp_type) THEN
            NULLIFY(geometric_parameters)
            CALL field_parameter_set_data_get(component%FIELD_VARIABLE%FIELD,component%FIELD_VARIABLE%VARIABLE_TYPE, &
              & field_values_set_type,geometric_parameters,err,error,*999)
            err = fieldexport_elementgridvalues( sessionhandle, isfirstvalueset, 1, &
              & geometric_parameters(component%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP))
          ELSE IF(component%FIELD_VARIABLE%DATA_TYPE==field_intg_type) THEN
            NULLIFY(geometric_parameters_intg)
            CALL field_parameter_set_data_get(component%FIELD_VARIABLE%FIELD,component%FIELD_VARIABLE%VARIABLE_TYPE, &
              & field_values_set_type,geometric_parameters_intg,err,error,*999)
            ALLOCATE(geometric_parameters_dp(SIZE(geometric_parameters_intg)))
            IF(err/=0) CALL flagerror("Could not allocate geometric parameters dp", err, error,*999 )
            geometric_parameters_dp(1:SIZE(geometric_parameters_intg))= &
              & REAL(geometric_parameters_intg(1:size(geometric_parameters_intg)))
            err = fieldexport_elementgridvalues( sessionhandle, isfirstvalueset, 1, &
              & geometric_parameters_dp(component%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP))
            DEALLOCATE(geometric_parameters_dp)
          ELSE
            CALL flagerror( "Only INTG and REAL data types implemented.", err, error,*999 )
          ENDIF
          isfirstvalueset = 0
        ELSE IF( component%INTERPOLATION_TYPE == field_gauss_point_based_interpolation) THEN
          IF(component%FIELD_VARIABLE%DATA_TYPE==field_dp_type) THEN
            NULLIFY(geometric_parameters)
            CALL field_parameter_set_data_get(component%FIELD_VARIABLE%FIELD,component%FIELD_VARIABLE%VARIABLE_TYPE, &
              & field_values_set_type,geometric_parameters,err,error,*999)
            err = fieldexport_elementgridvalues( sessionhandle, isfirstvalueset, basis%QUADRATURE%QUADRATURE_SCHEME_MAP( &
              & basis_default_quadrature_scheme)%PTR%NUMBER_OF_GAUSS, &
              & geometric_parameters(component%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(1,local_number)))
          ELSE IF(component%FIELD_VARIABLE%DATA_TYPE==field_intg_type) THEN
            NULLIFY(geometric_parameters_intg)
            CALL field_parameter_set_data_get(component%FIELD_VARIABLE%FIELD,component%FIELD_VARIABLE%VARIABLE_TYPE, &
              & field_values_set_type,geometric_parameters_intg,err,error,*999)
            ALLOCATE(geometric_parameters_dp(SIZE(geometric_parameters_intg)))
            IF(err/=0) CALL flagerror("Could not allocate geometric parameters dp", err, error,*999 )
            geometric_parameters_dp(1:SIZE(geometric_parameters_intg))= &
              & REAL(geometric_parameters_intg(1:size(geometric_parameters_intg)))
            err = fieldexport_elementgridvalues( sessionhandle, isfirstvalueset, basis%QUADRATURE%QUADRATURE_SCHEME_MAP( &
              & basis_default_quadrature_scheme)%PTR%NUMBER_OF_GAUSS, &
              & geometric_parameters_dp(component%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(1,local_number)))
            DEALLOCATE(geometric_parameters_dp)
          ELSE
            CALL flagerror( "Only INTG and REAL data types implemented.", err, error,*999 )
          ENDIF
          isfirstvalueset = 0
        ENDIF

        IF(err/=0) THEN
          CALL flagerror( "Cannot write grid points to nodes file", err, error,*999 )
        ENDIF
      ENDDO


      basis=>element%BASIS

!!TEMP
      !ERR = FieldExport_ElementNodeIndices( sessionHandle, BASIS%NUMBER_OF_NODES, C_LOC( element%USER_ELEMENT_NODES ) )
!!Copy user element nodes to a temporary array that has the target attribute. gcc bug 38813 prevents using C_LOC with
!!the array directly. nb using a fixed length array here which is dangerous but should suffice for now.
      SELECT CASE(basis%TYPE)
      CASE(basis_lagrange_hermite_tp_type)
        user_element_nodes(1:basis%NUMBER_OF_NODES)=element%USER_ELEMENT_NODES(1:basis%NUMBER_OF_NODES)
      CASE(basis_simplex_type)
        SELECT CASE(basis%NUMBER_OF_XI)
        CASE(1)
          user_element_nodes(1:basis%NUMBER_OF_NODES)=element%USER_ELEMENT_NODES(1:basis%NUMBER_OF_NODES)
        CASE(2)
          SELECT CASE(basis%INTERPOLATION_ORDER(1))
          CASE(basis_linear_interpolation_order)
            user_element_nodes(1:3)=element%USER_ELEMENT_NODES(1:3)
          CASE(basis_quadratic_interpolation_order)
            user_element_nodes(1)=element%USER_ELEMENT_NODES(1)
            user_element_nodes(2)=element%USER_ELEMENT_NODES(4)
            user_element_nodes(3)=element%USER_ELEMENT_NODES(2)
            user_element_nodes(4)=element%USER_ELEMENT_NODES(6)
            user_element_nodes(5)=element%USER_ELEMENT_NODES(5)
            user_element_nodes(6)=element%USER_ELEMENT_NODES(3)
          CASE(basis_cubic_interpolation_order)
            user_element_nodes(1)=element%USER_ELEMENT_NODES(1)
            user_element_nodes(2)=element%USER_ELEMENT_NODES(4)
            user_element_nodes(3)=element%USER_ELEMENT_NODES(5)
            user_element_nodes(4)=element%USER_ELEMENT_NODES(2)
            user_element_nodes(5)=element%USER_ELEMENT_NODES(9)
            user_element_nodes(6)=element%USER_ELEMENT_NODES(10)
            user_element_nodes(7)=element%USER_ELEMENT_NODES(6)
            user_element_nodes(8)=element%USER_ELEMENT_NODES(8)
            user_element_nodes(9)=element%USER_ELEMENT_NODES(7)
            user_element_nodes(10)=element%USER_ELEMENT_NODES(3)
         CASE DEFAULT
            CALL flagerror("Invalid basis order.",err,error,*999)
          END SELECT
        CASE(3)
          SELECT CASE(basis%INTERPOLATION_ORDER(1))
          CASE(basis_linear_interpolation_order)
            user_element_nodes(1:4)=element%USER_ELEMENT_NODES(1:4)
          CASE(basis_quadratic_interpolation_order)
            user_element_nodes(1)=element%USER_ELEMENT_NODES(1)
            user_element_nodes(2)=element%USER_ELEMENT_NODES(5)
            user_element_nodes(3)=element%USER_ELEMENT_NODES(2)
            user_element_nodes(4)=element%USER_ELEMENT_NODES(6)
            user_element_nodes(5)=element%USER_ELEMENT_NODES(8)
            user_element_nodes(6)=element%USER_ELEMENT_NODES(3)
            user_element_nodes(7)=element%USER_ELEMENT_NODES(7)
            user_element_nodes(8)=element%USER_ELEMENT_NODES(10)
            user_element_nodes(9)=element%USER_ELEMENT_NODES(9)
            user_element_nodes(10)=element%USER_ELEMENT_NODES(4)
          CASE(basis_cubic_interpolation_order)
            user_element_nodes(1)=element%USER_ELEMENT_NODES(1)
            user_element_nodes(2)=element%USER_ELEMENT_NODES(5)
            user_element_nodes(3)=element%USER_ELEMENT_NODES(6)
            user_element_nodes(4)=element%USER_ELEMENT_NODES(2)
            user_element_nodes(5)=element%USER_ELEMENT_NODES(7)
            user_element_nodes(6)=element%USER_ELEMENT_NODES(17)
            user_element_nodes(7)=element%USER_ELEMENT_NODES(11)
            user_element_nodes(8)=element%USER_ELEMENT_NODES(8)
            user_element_nodes(9)=element%USER_ELEMENT_NODES(12)
            user_element_nodes(10)=element%USER_ELEMENT_NODES(3)
            user_element_nodes(11)=element%USER_ELEMENT_NODES(9)
            user_element_nodes(12)=element%USER_ELEMENT_NODES(18)
            user_element_nodes(13)=element%USER_ELEMENT_NODES(15)
            user_element_nodes(14)=element%USER_ELEMENT_NODES(19)
            user_element_nodes(15)=element%USER_ELEMENT_NODES(20)
            user_element_nodes(16)=element%USER_ELEMENT_NODES(13)
            user_element_nodes(17)=element%USER_ELEMENT_NODES(10)
            user_element_nodes(18)=element%USER_ELEMENT_NODES(16)
            user_element_nodes(19)=element%USER_ELEMENT_NODES(14)
            user_element_nodes(20)=element%USER_ELEMENT_NODES(4)
          CASE DEFAULT
            CALL flagerror("Invalid basis order.",err,error,*999)
          END SELECT
        CASE DEFAULT
          CALL flagerror("Invalid number of xi.",err,error,*999)
        END SELECT
      CASE DEFAULT
        CALL flagerror("Not implemented.",err,error,*999)
      END SELECT
      err = fieldexport_elementnodeindices( sessionhandle, basis%NUMBER_OF_NODES, c_loc( user_element_nodes ) )
      IF(err/=0) THEN
        CALL flagerror( "Cannot write node indices to file", err, error,*999 )
      ENDIF

      CALL field_io_export_element_scale_factors( sessionhandle, components, &
        & list_comp_scale, global_number, my_computational_node_number, err, error, *999 )

    ENDDO !elem_idx

    err = fieldexport_closesession( sessionhandle )
    IF(err/=0) THEN
      CALL flagerror( "Cannot close element export file", err, error,*999 )
    ENDIF
    sessionhandle = -1

    !release the temporary memory
    CALL checked_deallocate( scale_factors )
    CALL checked_deallocate( nodal_number )
    CALL checked_deallocate( list_comp_scale )

    exits("FIELD_IO_EXPORT_ELEMENTS_INTO_LOCAL_FILE")
    RETURN
999 errorsexits("FIELD_IO_EXPORT_ELEMENTS_INTO_LOCAL_FILE",err,error)
    RETURN 1
  END SUBROUTINE field_io_export_elements_into_local_file

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

  SUBROUTINE field_io_elemental_info_set_sort(ELEMENTAL_INFO_SET, my_computational_node_number, ERR,ERROR,*)
    !Argument variables
    TYPE(field_io_info_set), INTENT(INOUT) :: ELEMENTAL_INFO_SET
    INTEGER(INTG), INTENT(IN):: my_computational_node_number
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(domain_mapping_type), POINTER :: DOMAIN_MAPPING_ELEMENTS !The domain mapping to calculate nodal mappings
    TYPE(domain_elements_type), POINTER :: DOMAIN_ELEMENTS1, DOMAIN_ELEMENTS2! domain nodes
    INTEGER(INTG) :: global_number1, local_number1, global_number2, local_number2
    INTEGER(INTG) :: component_idx, nn1, nn2 ! nn, tmp2, tmp1!temporary variable
    INTEGER(INTG) :: node_idx, deriv_idx
    TYPE(field_io_component_info_set), POINTER :: tmpInfoSet
    LOGICAL :: SAME_ELEMENT_INFO

    !from now on, global numbering are used
    enters("FIELD_IO_ELEMENTAL_INFO_SET_SORT",err,error,*999)

    IF(.NOT.ALLOCATED(elemental_info_set%LIST_OF_GLOBAL_NUMBER)) THEN
      CALL flagerror("list of global numbering in the input data is invalid",err,error,*999)
    ENDIF
    IF(.NOT.ALLOCATED(elemental_info_set%COMPONENT_INFO_SET)) THEN
      CALL flagerror("nodal information set in the input data is invalid",err,error,*999)
    ENDIF


    !!get my own computianal node number--be careful the rank of process in the MPI pool
    !!is not necessarily equal to numbering of computional node, so use method COMPUTATIONAL_NODE_NUMBER_GET
    !!will be a secured way to get the number
    !my_computational_node_number=COMPUTATIONAL_NODE_NUMBER_GET(ERR,ERROR)
    !IF(ERR/=0) GOTO 999

    !group nodal information set according to its components, i.e. put all the nodes with the same components together
    !and change the global number in the LIST_OF_GLOBAL_NUMBER
    nn1=1
    DO WHILE(nn1<elemental_info_set%NUMBER_OF_ENTRIES)
      !global number of this node
      global_number1=elemental_info_set%LIST_OF_GLOBAL_NUMBER(nn1)
      DO nn2=nn1+1,elemental_info_set%NUMBER_OF_ENTRIES
        global_number2=elemental_info_set%LIST_OF_GLOBAL_NUMBER(nn2)
        IF(elemental_info_set%COMPONENT_INFO_SET(nn1)%PTR%NUMBER_OF_COMPONENTS== &
            & elemental_info_set%COMPONENT_INFO_SET(nn2)%PTR%NUMBER_OF_COMPONENTS) THEN
          same_element_info=.true.
          !we will check the component (type of component, partial derivative).
          DO component_idx=1,elemental_info_set%COMPONENT_INFO_SET(nn1)%PTR%NUMBER_OF_COMPONENTS
            !not safe, but it is fast
            !=============================================================================================!
            !        checking according to local memory adddress                          !
            !=============================================================================================!
            !are they in the same memory address?
            IF(.NOT.ASSOCIATED(elemental_info_set%COMPONENT_INFO_SET(nn1)%PTR%COMPONENTS(component_idx)%PTR, &
              &TARGET=elemental_info_set%COMPONENT_INFO_SET(nn2)%PTR%COMPONENTS(component_idx)%PTR))  THEN
              same_element_info=.false.
              EXIT
            ENDIF !ASSCOCIATED

            !! better use this one because it is safe method, but slow
            !!=============================================================================================!
            !!        checking according to the types defined in the openCMISS                  !
            !!=============================================================================================!
            !!are they in the same field?
            !IF(ELEMENTAL_INFO_SET%COMPONENT_INFO_SET(nn1)%PTR%COMPONENTS(component_idx)%PTR%FIELD%GLOBAL_NUMBER/= &
            !&ELEMENTAL_INFO_SET%COMPONENT_INFO_SET(nn2)%PTR%COMPONENTS(component_idx)%PTR%FIELD%GLOBAL_NUMBER) THEN
            !  SAME_ELEMENT_INFO=.FALSE.
            !  EXIT
            !ELSE  !GLOBAL_NUBMER
            !  !are they the same variable?
            !  IF(ELEMENTAL_INFO_SET%COMPONENT_INFO_SET(nn1)%PTR%COMPONENTS(component_idx)%PTR%FIELD_VARIABLE%VARIABLE_NUMBER/= &
            !  & ELEMENTAL_INFO_SET%COMPONENT_INFO_SET(nn2)%PTR%COMPONENTS(component_idx)%PTR%FIELD_VARIABLE%VARIABLE_NUMBER) THEN
            !     SAME_ELEMENT_INFO=.FALSE.
            !     EXIT
            !   ELSE !VARIABLE_NUBMER
            !    !are they the same component?
            !    IF(LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn1)%PTR%COMPONENTS(component_idx)%PTR%COMPONENT_NUMBER/=&
            !      &LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn2)%PTR%COMPONENTS(component_idx)%PTR%COMPONENT_NUMBER) THEN
            !       SAME_ELEMENT_INFO=.FALSE.
            !       EXIT
            !     ENDIF !COMPONENT_NUMBER
            !  ENDIF ! VARIABLE_NUBMER
            !ENDIF !GLOBAL_NUBMER
          ENDDO !component_idx

          !check whether correspoding two components have the same partial derivatives
          IF(same_element_info) THEN
            DO component_idx=1,elemental_info_set%COMPONENT_INFO_SET(nn1)%PTR%NUMBER_OF_COMPONENTS
              !finding the local numbering for the NODAL_INFO_SET(nn1)
              domain_mapping_elements=>&
                & elemental_info_set%COMPONENT_INFO_SET(nn1)%PTR%COMPONENTS(component_idx)%PTR%DOMAIN% &
                & mappings%ELEMENTS
              !get the domain index for this variable component according to my own computional node number
              !local number of nn1'th node in the damain assoicated with component(component_idx)
              local_number1 = findmylocaldomainnumber( domain_mapping_elements%GLOBAL_TO_LOCAL_MAP( global_number1 ), &
                & my_computational_node_number )
              domain_elements1=>&
                & elemental_info_set%COMPONENT_INFO_SET(nn1)%PTR%COMPONENTS(component_idx)%PTR% &
                & domain%TOPOLOGY%ELEMENTS

              !finding the local numbering for the NODAL_INFO_SET(nn2)
              domain_mapping_elements=>&
                & elemental_info_set%COMPONENT_INFO_SET(nn2)%PTR%COMPONENTS(component_idx)%PTR% &
                & domain%MAPPINGS%ELEMENTS
              !get the domain index for this variable component according to my own computional node number
              !local number of nn2'th node in the damain assoicated with component(component_idx)
              local_number2 = findmylocaldomainnumber( domain_mapping_elements%GLOBAL_TO_LOCAL_MAP( global_number2 ), &
                & my_computational_node_number )
              domain_elements2=>&
                & elemental_info_set%COMPONENT_INFO_SET(nn2)%PTR%COMPONENTS(component_idx)%PTR% &
                & domain%TOPOLOGY%ELEMENTS

              !checking whether they have the same basis
              IF(domain_elements1%ELEMENTS(local_number1)%BASIS%GLOBAL_NUMBER/=&
                &domain_elements2%ELEMENTS(local_number2)%BASIS%GLOBAL_NUMBER) THEN
                same_element_info=.false.
                EXIT
              ENDIF  !DOMAIN_ELEMENTS1

              !Check that the elements use the same versions of all derivatives for all element nodes
              DO node_idx=1,domain_elements1%ELEMENTS(local_number1)%BASIS%NUMBER_OF_NODES
                DO deriv_idx=1,domain_elements1%ELEMENTS(local_number1)%BASIS%NUMBER_OF_DERIVATIVES(node_idx)
                  IF (domain_elements1%ELEMENTS(local_number1)%elementVersions(deriv_idx,node_idx)/= &
                      & domain_elements2%ELEMENTS(local_number2)%elementVersions(deriv_idx,node_idx)) THEN
                    same_element_info=.false.
                    EXIT
                  END IF
                END DO
              END DO

            ENDDO !component_idx
          ENDIF !SAME_ELEMENT_INFO==.TRUE.
        ENDIF !LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS==LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn+1)%PTR%NUMBER_OF_COMPONENTS

        !find two elements which have the same output, and then they should put together
        IF(same_element_info) THEN
          tmpinfoset => elemental_info_set%COMPONENT_INFO_SET(nn2)%PTR
          elemental_info_set%COMPONENT_INFO_SET(nn2)%PTR => elemental_info_set%COMPONENT_INFO_SET(nn1+1)%PTR
          elemental_info_set%COMPONENT_INFO_SET(nn1+1)%PTR => tmpinfoset

          elemental_info_set%COMPONENT_INFO_SET(nn2)%PTR%SAME_HEADER=.false.
          elemental_info_set%COMPONENT_INFO_SET(nn1+1)%PTR%SAME_HEADER=.true.

          !exchange the global number
          elemental_info_set%LIST_OF_GLOBAL_NUMBER(nn2)=elemental_info_set%LIST_OF_GLOBAL_NUMBER(nn1+1)
          elemental_info_set%LIST_OF_GLOBAL_NUMBER(nn1+1)=global_number2

          !increase nn1 to skip the nodes which have the same output
          nn1=nn1+1
        ENDIF !(SAME_ELEMENT_INFO=.TRUE.)
      ENDDO !nn2
      !increase the nn1 to check next node
      nn1=nn1+1
    ENDDO !nn1<LOCAL_PROCESS_NODAL_INFO_SET%NUMBER_OF_ENTRIES

    !order the variable components and group them: X1(1),X1(2),X1(3),X2(2),X2(3),X3(2)....
    !DO nn=1,LOCAL_PROCESS_NODAL_INFO_SET%NUMBER_OF_ENTRIES
    !   print "(A, I)", "nn=", nn
    !   !temporarily use nk, nu here to save memory
    !   IF(LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS/=1) THEN
    !     component_idx=1
    !     DO WHILE(component_idx<LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS)
    !        !checking the same variable's components
    !       print "(A, I)", "component_idx=", component_idx
    !       print "(A, I)", "LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS", LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS
    !       DO WHILE(ASSOCIATED(LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS(component_idx)%PTR%FIELD_VARIABLE, &
    !       & TARGET=LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS(component_idx+1)%PTR%FIELD_VARIABLE))
    !          component_idx=component_idx+1
    !          IF(component_idx>=LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS) THEN
    !             EXIT
    !          ENDIF
    !        ENDDO
    !
    !        !It may have more than 3 component in the future?!! I do not know,too
    !        !so there the components are sorted according their numbering of component
    !        !nk and nu are used here temporarily
    !        DO tmp1=1,component_idx
    !           print "(A, I)", "tmp1=", tmp1
    !           SAME_ELEMENT_INFO=.FALSE.
    !           DO tmp2=1,(component_idx-tmp1)
    !              IF(LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS(tmp2)%PTR%COMPONENT_NUMBER>&
    !              &LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS(tmp2+1)%PTR%COMPONENT_NUMBER) THEN
    !                 tmp_ptr=>LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS(tmp2+1)%PTR
    !                 LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS(tmp2+1)%PTR=>&
    !                 LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS(tmp2)%PTR
    !
    !                 LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS(tmp2)%PTR=>tmp_ptr
    !                 SAME_ELEMENT_INFO=.TRUE.
    !              ENDIF
    !           ENDDO
    !           IF(SAME_ELEMENT_INFO) THEN
    !              EXIT
    !           ENDIF
    !        ENDDO
    !        NULLIFY(tmp_ptr)
    !        component_idx=component_idx+1
    !     ENDDO ! WHILE(component_idx<LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS)
    !   ENDIF ! LOCAL_PROCESS_NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS/=1
    !ENDDO !nn

    exits("FIELD_IO_ELEMENTAL_INFO_SET_SORT")
    RETURN
999 errorsexits("FIELD_IO_ELEMENTAL_INFO_SET_SORT",err,error)
    RETURN 1
  END SUBROUTINE field_io_elemental_info_set_sort

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

  SUBROUTINE fieldio_elementalinfosetattachlocalprocess( ELEMENTAL_INFO_SET, FIELDS, ERR, ERROR, * )
    !Argument variables
    TYPE(field_io_info_set), INTENT(INOUT):: ELEMENTAL_INFO_SET
    TYPE(fields_type), POINTER ::FIELDS
    INTEGER(INTG), INTENT(OUT):: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    LOGICAL :: ININTERFACE,INREGION
    TYPE(varying_string) :: LOCAL_ERROR
    TYPE(field_type), POINTER :: FIELD
    TYPE(domain_mapping_type), POINTER:: DOMAIN_ELEMENTS_MAPPING !nodes in local mapping--it is different as exnode
    TYPE(field_variable_type), POINTER:: FIELD_VARIABLE !field variable
    INTEGER(INTG) :: num_field, var_idx, component_idx, np, nn !temporary variable
    LOGICAL :: foundNewElement

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

    !validate the input data
    inregion=.false.
    ininterface=.false.
    IF(ASSOCIATED(fields%REGION)) THEN
      inregion=.true.
    ELSE
      IF(ASSOCIATED(fields%INTERFACE)) THEN
        ininterface=.true.
      ELSE
        CALL flagerror("Fields is not associated with a region or interface.",err,error,*999)
      ENDIF
    ENDIF

    IF(inregion) THEN
      !checking whether the list of fields in the same region
      DO num_field =1, fields%NUMBER_OF_FIELDS
        IF(.NOT.ASSOCIATED(fields%FIELDS(num_field)%PTR)) THEN
          local_error ="No. "//trim(number_to_vstring(num_field,"*",err,error))// &
            & " field handle in fields list is invalid"
          CALL flagerror(local_error,err,error,*999)
        ENDIF

        IF( num_field == 1 ) THEN
          cycle
        ENDIF

        IF(fields%FIELDS(num_field-1)%PTR%REGION%USER_NUMBER/=fields%FIELDS(num_field)%PTR%REGION%USER_NUMBER) THEN
          local_error = "No. "//trim(number_to_vstring(num_field-1,"*",err,error))//" and "// &
            & trim(number_to_vstring(num_field,"*",err,error))//" fields are not in the same region"
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ENDDO
    ELSE
      !checking whether the list of fields in the same interface
      DO num_field =1, fields%NUMBER_OF_FIELDS
        IF(.NOT.ASSOCIATED(fields%FIELDS(num_field)%PTR)) THEN
          local_error ="No. "//trim(number_to_vstring(num_field,"*",err,error))// &
            & " field handle in fields list is invalid"
          CALL flagerror(local_error,err,error,*999)
        ENDIF

        IF( num_field == 1 ) THEN
          cycle
        ENDIF

        IF(fields%FIELDS(num_field-1)%PTR%INTERFACE%USER_NUMBER/= &
          & fields%FIELDS(num_field)%PTR%INTERFACE%USER_NUMBER) THEN
          local_error = "No. "//trim(number_to_vstring(num_field-1,"*",err,error))//" and "// &
            & trim(number_to_vstring(num_field,"*",err,error))//" fields are not in the same interface."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ENDDO
    ENDIF

    elemental_info_set%FIELDS=>fields

    !attache local process to local nodal information set. In current opencmiss system,
    !each local process owns it local nodal information, so all we need to do is to fill the nodal
    !information set with nodal information of local process
    IF((elemental_info_set%NUMBER_OF_ENTRIES/=0).OR.(.NOT.ASSOCIATED(elemental_info_set%FIELDS)) &
      & .OR.ALLOCATED(elemental_info_set%COMPONENT_INFO_SET)) THEN
      CALL flagerror("nodal information set is not initialized properly, and call start method first", &
        & err,error,*999)
    ENDIF

    DO num_field=1,elemental_info_set%FIELDS%NUMBER_OF_FIELDS
      field=>elemental_info_set%FIELDS%FIELDS(num_field)%PTR
      IF(.NOT.ALLOCATED(field%VARIABLES)) THEN
        cycle
      ENDIF
      DO var_idx=1, field%NUMBER_OF_VARIABLES
        field_variable=>field%VARIABLES(var_idx)
        DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
          IF(.NOT.ASSOCIATED(field_variable%COMPONENTS(component_idx)%DOMAIN%TOPOLOGY%ELEMENTS)) THEN
            cycle
          ENDIF

          domain_elements_mapping=>field_variable%COMPONENTS(component_idx)%DOMAIN%MAPPINGS%ELEMENTS
          DO np=1,domain_elements_mapping%NUMBER_OF_LOCAL
            foundnewelement=.true.
            DO nn=1,elemental_info_set%NUMBER_OF_ENTRIES
              IF(elemental_info_set%LIST_OF_GLOBAL_NUMBER(nn)==domain_elements_mapping%LOCAL_TO_GLOBAL_MAP(np)) THEN
                foundnewelement=.false.
                EXIT
              ENDIF
            ENDDO
            !have one more global node
            !i hate the codes here, but i have to save the memory
            IF(foundnewelement) THEN
              CALL grow_array( elemental_info_set%LIST_OF_GLOBAL_NUMBER, 1, &
                & "Could not allocate temporary buffer in IO", err, error, *999 )
              elemental_info_set%LIST_OF_GLOBAL_NUMBER(elemental_info_set% &
                & number_of_entries+1) = domain_elements_mapping%LOCAL_TO_GLOBAL_MAP(np)
              elemental_info_set%NUMBER_OF_ENTRIES=elemental_info_set%NUMBER_OF_ENTRIES+1
            ENDIF !foundNewElement
          ENDDO !np
        ENDDO !component_idx
      ENDDO !var_idx
    ENDDO !num_field

    !allocate the nodal information set and initialize them
    ALLOCATE(elemental_info_set%COMPONENT_INFO_SET(elemental_info_set%NUMBER_OF_ENTRIES),stat=err)
    IF(err/=0) CALL flagerror("Could not allocate nodal information set",err,error,*999)

    DO nn = 1, elemental_info_set%NUMBER_OF_ENTRIES
      ALLOCATE( elemental_info_set%COMPONENT_INFO_SET(nn)%PTR )
      elemental_info_set%COMPONENT_INFO_SET(nn)%PTR%SAME_HEADER = .false.
      elemental_info_set%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS = 0
      CALL checked_deallocate( elemental_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS )
    ENDDO

    !collect nodal information from local process
    DO num_field=1,elemental_info_set%FIELDS%NUMBER_OF_FIELDS
      field=>elemental_info_set%FIELDS%FIELDS(num_field)%PTR
      IF(.NOT.ALLOCATED(field%VARIABLES)) THEN
        cycle
      ENDIF
      DO var_idx=1, field%NUMBER_OF_VARIABLES
        field_variable=>field%VARIABLES(var_idx)
        DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
          IF(.NOT.ASSOCIATED(field_variable%COMPONENTS(component_idx)%DOMAIN%TOPOLOGY%ELEMENTS)) THEN
            cycle
          ENDIF

          domain_elements_mapping=>field_variable%COMPONENTS(component_idx)%DOMAIN%MAPPINGS%ELEMENTS
          DO np=1,domain_elements_mapping%NUMBER_OF_LOCAL
            DO nn=1,elemental_info_set%NUMBER_OF_ENTRIES
              IF(elemental_info_set%LIST_OF_GLOBAL_NUMBER(nn)==domain_elements_mapping%LOCAL_TO_GLOBAL_MAP(np)) THEN
                EXIT
              ENDIF
            ENDDO

            !allocate variable component memory
            CALL grow_array( elemental_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS, 1, &
              & "Could not allocate component buffer in IO", err, error, *999 )
            elemental_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS( &
              & elemental_info_set%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS+1 &
              & )%PTR=>field%VARIABLES(var_idx)%COMPONENTS(component_idx)
            !increase number of component
            elemental_info_set%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS=&
              & elemental_info_set%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS+1
          ENDDO !np
        ENDDO !component_idx
      ENDDO !var_idx
    ENDDO !num_field

    !LOCAL_PROCESS_NODAL_INFO_SET%LIST_OF_GLOBAL_NUMBER=>LIST_OF_GLOBAL_NUMBER
    !NULLIFY(LIST_OF_GLOBAL_NUMBER)

    exits("FieldIO_ElementalInfoSetAttachLocalProcess")
    RETURN
999 errorsexits("FieldIO_ElementalInfoSetAttachLocalProcess",err,error)
    RETURN 1
  END SUBROUTINE fieldio_elementalinfosetattachlocalprocess

  !
  !================================================================================================================================
  !
  !
  !!>Import nodal information \see{FIELD_IO::FIELD_IO_NODES_IMPORT}.
  !SUBROUTINE FIELD_IO_NODES_IMPORT(FIELDS, FILE_NAME, METHOD, ERR,ERROR,*)
  !  !Argument variables
  !  TYPE(FIELDS_TYPE), POINTER :: FIELDS !<the field object
  !  TYPE(VARYING_STRING), INTENT(INOUT) :: FILE_NAME !<file name
  !  TYPE(VARYING_STRING), INTENT(IN):: METHOD
  !  INTEGER(INTG), INTENT(OUT) :: ERR !<The error code
  !  TYPE(VARYING_STRING), INTENT(OUT) :: ERROR !<The error string
  !  !Local Variables
  !  TYPE(FIELD_IO_INFO_SET) :: LOCAL_PROCESS_NODAL_INFO_SET !<nodal information in this process
  !  INTEGER(INTG):: my_computational_node_number !<local process number
  !  INTEGER(INTG):: computational_node_numbers   !<total process number
  !
  !  ENTERS("FIELD_IO_NODES_IMPORT", ERR,ERROR,*999)
  !
  !  !Get the number of computational nodes
  !  computational_node_numbers=COMPUTATIONAL_NODES_NUMBER_GET(ERR,ERROR)
  !  IF(ERR/=0) GOTO 999
  !  !Get my computational node number
  !  my_computational_node_number=COMPUTATIONAL_NODE_NUMBER_GET(ERR,ERROR)
  !  IF(ERR/=0) GOTO 999
  !  IF(METHOD=="FORTRAN") THEN
  !     CALL FIELD_IO_INFO_SET_INITIALISE(LOCAL_PROCESS_NODAL_INFO_SET, FIELDS, ERR,ERROR,*999)
  !     CALL FieldIO_NodelInfoSetAttachLocalProcess(LOCAL_PROCESS_NODAL_INFO_SET, ERR,ERROR,*999)
  !     CALL FIELD_IO_NODAL_INFO_SET_SORT(LOCAL_PROCESS_NODAL_INFO_SET, my_computational_node_number, ERR,ERROR,*999)
  !     CALL FIELD_IO_IMPORT_NODES_FROM_LOCAL_FILE(LOCAL_PROCESS_NODAL_INFO_SET, FILE_NAME, my_computational_node_number, &
  !          &computational_node_numbers, ERR, ERROR, *999)
  !     CALL FIELD_IO_NODAL_INFO_SET_FINALIZE(LOCAL_PROCESS_NODAL_INFO_SET, ERR,ERROR,*999)
  !  ELSE IF(METHOD=="MPIIO") THEN
  !     CALL FlagError("what are u thinking, of course not!",ERR,ERROR,*999)
  !  ENDIF
  !
  !  EXITS("FIELD_IO_NODES_IMPORT")
  !  RETURN
!999 ERRORSEXITS("FIELD_IO_NODES_IMPORT",ERR,ERROR)
  !  RETURN 1
  !END SUBROUTINE FIELD_IO_NODES_IMPORT


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

  FUNCTION field_io_compare_info_set_components( SET1, SET2 )
    !Argument variables
    TYPE(field_io_component_info_set) :: SET1
    TYPE(field_io_component_info_set) :: SET2
    LOGICAL :: FIELD_IO_COMPARE_INFO_SET_COMPONENTS

    !local variables
    INTEGER(INTG) :: component_idx

    field_io_compare_info_set_components = .false.

    IF( set1%NUMBER_OF_COMPONENTS /= set2%NUMBER_OF_COMPONENTS ) THEN
      RETURN
    ENDIF

    DO component_idx = 1, set1%NUMBER_OF_COMPONENTS
      !!not safe, but it is fast
      !!=============================================================================================!
      !!           checking according to local memory adddress                                       !
      !!=============================================================================================!
      !!are they in the same memory address?
      !IF(SET1%COMPONENTS(component_idx)%PTR/=&
      !  &SET2%COMPONENTS(component_idx)%PTR)
      !THEN
      !   FIELD_IO_COMPARE_INFO_SETS=.FALSE. !out of loop-component_idx=1,SET1%NUMBER_OF_COMPONENTS
      !   EXIT
      !ENDIF                      NUMBER_OF_NODES

      ! better use this one because it is safe method, but slow
      !=============================================================================================!
      !           checking according to the types defined in the openCMISS                          !
      !=============================================================================================!
      !are they in the same field?
      IF( set1%COMPONENTS(component_idx)%PTR%FIELD_VARIABLE%FIELD%GLOBAL_NUMBER/= &
        & set2%COMPONENTS(component_idx)%PTR%FIELD_VARIABLE%FIELD%GLOBAL_NUMBER ) THEN
        RETURN
      ENDIF

        !are they the same variable?
      IF( set1%COMPONENTS( component_idx )%PTR%FIELD_VARIABLE% &
        & variable_number /= set2%COMPONENTS( component_idx )%PTR% &
        & field_variable%VARIABLE_NUMBER ) THEN
        RETURN
      ENDIF

      !are they the same component?
      IF( set1%COMPONENTS( component_idx )%PTR%COMPONENT_NUMBER /= &
        & set2%COMPONENTS( component_idx)%PTR%COMPONENT_NUMBER ) THEN
        RETURN
      ENDIF
    ENDDO !component_idx

    field_io_compare_info_set_components = .true.

  END FUNCTION field_io_compare_info_set_components

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

  SUBROUTINE field_io_compare_info_set_derivatives( SET1, SET2, my_computational_node_number, global_number1, global_number2, &
    & doesmatch, err, error, * )
    !Argument variables
    TYPE(field_io_component_info_set) :: SET1
    TYPE(field_io_component_info_set) :: SET2
    INTEGER(INTG) :: my_computational_node_number
    INTEGER(INTG) :: global_number1
    INTEGER(INTG) :: global_number2
    LOGICAL :: doesMatch
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    !local variables
    INTEGER(INTG) :: component_idx, derivative_idx
    INTEGER(INTG) :: local_number1, local_number2, tmp1
    TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES1, DOMAIN_NODES2
    INTEGER(INTG), ALLOCATABLE:: array1(:), array2(:)
    LOGICAL :: FOUND

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

    doesmatch = .true.

    !We have a potential match. Do a deeper inspection
    DO component_idx=1, set1%NUMBER_OF_COMPONENTS

      domain_nodes1=>set1%COMPONENTS(component_idx)%PTR%DOMAIN%TOPOLOGY%NODES
      found=.false.
      DO local_number1=1,domain_nodes1%NUMBER_OF_NODES
        IF( domain_nodes1%NODES(local_number1)%GLOBAL_NUMBER == global_number1 ) THEN
          found = .true.
          EXIT
        ENDIF
      ENDDO !local_number

      IF( .NOT. found ) THEN
        doesmatch = .false.
        EXIT !out of loop-component_idx=1,SET1%NUMBER_OF_COMPONENTS
      ENDIF

      domain_nodes2=>set2%COMPONENTS(component_idx)%PTR%DOMAIN%TOPOLOGY%NODES
      found=.false.
      DO local_number2=1,domain_nodes2%NUMBER_OF_NODES
        IF( domain_nodes2%NODES(local_number2)%GLOBAL_NUMBER == global_number2 ) THEN
          found = .true.
          EXIT
        ENDIF
      ENDDO !local_number

      IF( .NOT. found ) THEN
        doesmatch = .false.
        EXIT !out of loop-component_idx=1,SET1%NUMBER_OF_COMPONENTS
      ENDIF

      IF(domain_nodes1%NODES(local_number1)%NUMBER_OF_DERIVATIVES&
        &==domain_nodes2%NODES(local_number2)%NUMBER_OF_DERIVATIVES) THEN
        ALLOCATE(array1(domain_nodes1%NODES(local_number1)%NUMBER_OF_DERIVATIVES),stat=err)
        IF(err/=0) CALL flagerror("Could not allocate temporary buffer in IO sorting",err,error,*999)

        ALLOCATE(array2(domain_nodes1%NODES(local_number2)%NUMBER_OF_DERIVATIVES),stat=err)
        IF(err/=0) CALL flagerror("Could not allocate temporary buffer in IO sorting",err,error,*999)

        array1(1:domain_nodes1%NODES(local_number1)%NUMBER_OF_DERIVATIVES)=0
        array2(1:domain_nodes1%NODES(local_number2)%NUMBER_OF_DERIVATIVES)=0

        DO derivative_idx=1,domain_nodes1%NODES(local_number1)%NUMBER_OF_DERIVATIVES
          array1(derivative_idx)=domain_nodes1%NODES(local_number1)%DERIVATIVES(derivative_idx)%PARTIAL_DERIVATIVE_INDEX
        ENDDO
        DO derivative_idx=1,domain_nodes1%NODES(local_number2)%NUMBER_OF_DERIVATIVES
          array2(derivative_idx)=domain_nodes1%NODES(local_number2)%DERIVATIVES(derivative_idx)%PARTIAL_DERIVATIVE_INDEX
        ENDDO

        CALL list_sort(array1,err,error,*999)
        CALL list_sort(array2,err,error,*999)
        tmp1=sum(array1-array2)
        DEALLOCATE(array1)
        DEALLOCATE(array2)
        IF(tmp1/=0) THEN
          doesmatch = .false.
          EXIT !out of loop-component_idx=1,SET1%NUMBER_OF_COMPONENTS
        ENDIF
      ELSE
        doesmatch = .false.
        EXIT
      ENDIF

      ! Check that the nodes have the same number of versions, otherwise they must be grouped separately
      IF(set1%COMPONENT_VERSIONS(component_idx)/=set2%COMPONENT_VERSIONS(component_idx)) THEN
        doesmatch = .false.
        EXIT
      END IF
    ENDDO !component_idx

    exits("FIELD_IO_COMPARE_INFO_SET_DERIVATIVES")
    RETURN
999 errorsexits("FIELD_IO_COMPARE_INFO_SET_DERIVATIVES",err,error)
    RETURN 1

  END SUBROUTINE field_io_compare_info_set_derivatives

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

  SUBROUTINE field_io_nodal_info_set_sort(NODAL_INFO_SET, my_computational_node_number, ERR,ERROR,*)
    !Argument variables
    TYPE(field_io_info_set), INTENT(INOUT) :: NODAL_INFO_SET
    INTEGER(INTG), INTENT(IN):: my_computational_node_number
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(field_io_component_info_set), POINTER :: tmpInfoSet
    INTEGER(INTG) :: global_number1, global_number2
    INTEGER(INTG) :: nn1, nn2
    LOGICAL :: SAME_NODAL_INFO

    !from now on, global numbering are used
    enters("FIELD_IO_NODAL_INFO_SET_SORT",err,error,*999)

    IF(.NOT.ALLOCATED(nodal_info_set%LIST_OF_GLOBAL_NUMBER)) THEN
       CALL flagerror("list of global numbering in the input data is invalid",err,error,*999)
    ENDIF
    IF(.NOT.ALLOCATED(nodal_info_set%COMPONENT_INFO_SET)) THEN
       CALL flagerror("nodal information set in the input data is invalid",err,error,*999)
    ENDIF

    !group nodal information set according to its components, i.e. put all the nodes with the same components together
    !and change the global number in the LIST_OF_GLOBAL_NUMBER
    nn1=1
    DO WHILE(nn1<nodal_info_set%NUMBER_OF_ENTRIES)
      !global number of this node
      global_number1=nodal_info_set%LIST_OF_GLOBAL_NUMBER(nn1)
      DO nn2=nn1+1,nodal_info_set%NUMBER_OF_ENTRIES
        global_number2=nodal_info_set%LIST_OF_GLOBAL_NUMBER(nn2)

        same_nodal_info = field_io_compare_info_set_components( nodal_info_set%COMPONENT_INFO_SET( nn1 )%PTR, &
          & nodal_info_set%COMPONENT_INFO_SET( nn2 )%PTR )

        !check whether correspoding two components have the same partial derivatives
        IF( same_nodal_info ) THEN
          CALL field_io_compare_info_set_derivatives( nodal_info_set%COMPONENT_INFO_SET(nn1)%PTR, &
              & nodal_info_set%COMPONENT_INFO_SET(nn2)%PTR, my_computational_node_number, global_number1, global_number2, &
              & same_nodal_info, err, error, *999 )
        ENDIF !SAME_NODAL_INFO==.TRUE.

        !find two nodes which have the same output, and then they should put together
        IF(same_nodal_info) THEN
          tmpinfoset => nodal_info_set%COMPONENT_INFO_SET(nn2)%PTR
          nodal_info_set%COMPONENT_INFO_SET(nn2)%PTR => nodal_info_set%COMPONENT_INFO_SET(nn1+1)%PTR
          nodal_info_set%COMPONENT_INFO_SET(nn1+1)%PTR => tmpinfoset

          nodal_info_set%COMPONENT_INFO_SET(nn2)%PTR%SAME_HEADER=.false.
          nodal_info_set%COMPONENT_INFO_SET(nn1+1)%PTR%SAME_HEADER=.true.

          !exchange the global number
          nodal_info_set%LIST_OF_GLOBAL_NUMBER(nn2)=nodal_info_set%LIST_OF_GLOBAL_NUMBER(nn1+1)
          nodal_info_set%LIST_OF_GLOBAL_NUMBER(nn1+1)=global_number2

          !increase nn1 to skip the nodes which have the same output
          nn1=nn1+1
        ENDIF !(SAME_NODAL_INFO=.TRUE.)

      ENDDO !nn2
      !increase the nn1 to check next node
      nn1=nn1+1
    ENDDO !nn1<NODAL_INFO_SET%NUMBER_OF_ENTRIES

    !order the variable components and group them: X1(1),X1(2),X1(3),X2(2),X2(3),X3(2)....
    !DO nn=1,NODAL_INFO_SET%NUMBER_OF_ENTRIES
    !   print "(A, I)", "nn=", nn
    !   !temporarily use nk, nu here to save memory
    !   IF(NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%NUMBER_OF_COMPONENTS/=1) THEN
    !     component_idx=1
    !     DO WHILE(component_idx<NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%NUMBER_OF_COMPONNode:ENTS)
    !        !checking the same variable's components
    !       print "(A, I)", "component_idx=", component_idx
    !       print "(A, I)", "NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%NUMBER_OF_COMPONENTS", NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%NUMBER_OF_COMPONENTS
    !       DO WHILE(ASSOCIATED(NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%COMPONENTS(component_idx)%PTR%FIELD_VARIABLE, &
    !       & TARGET=NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%COMPONENTS(component_idx+1)%PTR%FIELD_VARIABLE))
    !          component_idx=component_idx+1
    !          IF(component_idx>=NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%NUMBER_OF_COMPONENTS) THEN
    !             EXIT
    !          ENDIF
    !        ENDDO
    !
    !        !It may have more than 3 component in the future?!! I do not know,too
    !        !so there the components are sorted according their numbering of component
    !        !nk and nu are used here temporarily
    !        DO tmp1=1,component_idx
    !           print "(A, I)", "tmp1=", tmp1
    !           SAME_NODAL_INFO=.FALSE.
    !           DO tmp2=1,(component_idx-tmp1)
    !              IF(NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%COMPONENTS(tmp2)%PTR%COMPONENT_NUMBER>&
    !              &NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%COMPONENTS(tmp2+1)%PTR%COMPONENT_NUMBER) THEN
    !                 tmp_ptr=>NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%COMPONENTS(tmp2+1)%PTR
    !                 NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%COMPONENTS(tmp2+1)%PTR=>&
    !                 NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%COMPONENTS(tmp2)%PTR
    !
    !                 NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%COMPONENTS(tmp2)%PTR=>tmp_ptr
    !                 SAME_NODAL_INFO=.TRUE.
    !              ENDIF
    !           ENDDO
    !           IF(SAME_NODAL_INFO) THEN
    !              EXIT
    !           ENDIF
    !        ENDDO
    !        NULLIFY(tmp_ptr)
    !        component_idx=component_idx+1
    !     ENDDO ! WHILE(component_idx<NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%NUMBER_OF_COMPONENTS)
    !   ENDIF ! NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%NUMBER_OF_COMPONENTS/=1
    !ENDDO !nn

    exits("FIELD_IO_NODAL_INFO_SET_SORT")
    RETURN
999 errorsexits("FIELD_IO_NODAL_INFO_SET_SORT",err,error)
    RETURN 1
  END SUBROUTINE field_io_nodal_info_set_sort

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

  FUNCTION field_io_label_derivative_info_get(GROUP_DERIVATIVES, NUMBER_DERIVATIVES, LABEL_TYPE, ERR, ERROR)
    !Argument variables
    INTEGER(INTG), INTENT(IN) :: NUMBER_DERIVATIVES
    INTEGER(INTG), INTENT(IN) :: GROUP_DERIVATIVES(number_derivatives)
    INTEGER(INTG), INTENT(IN) :: LABEL_TYPE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) ::FIELD_IO_LABEL_DERIVATIVE_INFO_GET
    INTEGER(INTG) :: dev_idx

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

    IF(number_derivatives==0) THEN
       CALL flagerror("number of derivatives in the input data is zero",err,error,*999)
    ENDIF
    IF(label_type/=field_io_derivative_label) THEN
       CALL flagerror("label type in the input data is not derivative label",err,error,*999)
    ENDIF

    IF((number_derivatives==1).AND.group_derivatives(1)==no_part_deriv) THEN
      field_io_label_derivative_info_get=" "
    ELSE
      field_io_label_derivative_info_get="("
      DO dev_idx=1,number_derivatives
        SELECT CASE(group_derivatives(dev_idx))
        CASE(no_part_deriv)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get
        CASE(part_deriv_s1)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d/ds1"
        CASE(part_deriv_s1_s1)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d2/ds1ds1"
        CASE(part_deriv_s2)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d/ds2"
        CASE(part_deriv_s2_s2)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d2/ds2ds2"
        CASE(part_deriv_s1_s2)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d/ds3"
        CASE(part_deriv_s3)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d2/ds3ds3"
        CASE(part_deriv_s3_s3)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d2/ds3ds3"
        CASE(part_deriv_s1_s3)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d2/ds1ds3"
        CASE(part_deriv_s2_s3)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d2/ds2ds3"
        CASE(part_deriv_s1_s2_s3)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d3/ds1ds2ds3"
        CASE(part_deriv_s4)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d/ds4"
        CASE(part_deriv_s4_s4)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d2/ds4ds4"
        CASE(part_deriv_s1_s4)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d2/ds1ds4"
        CASE(part_deriv_s2_s4)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d2/ds2ds4"
        CASE(part_deriv_s3_s4)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d2/ds3ds4"
        CASE(part_deriv_s1_s2_s4)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d3/ds1ds2ds4"
        CASE(part_deriv_s1_s3_s4)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d3/ds1ds3ds4"
        CASE(part_deriv_s2_s3_s4)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d3/ds2ds3ds4"
        CASE(part_deriv_s1_s4_s4)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d3/ds1ds4ds4"
        CASE(part_deriv_s2_s4_s4)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d3/ds2ds4ds4"
        CASE(part_deriv_s3_s4_s4)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d3/ds3ds4ds4"
        CASE(part_deriv_s4_s4_s4)
          field_io_label_derivative_info_get=field_io_label_derivative_info_get//", d3/ds4ds4ds4"
        CASE DEFAULT
          field_io_label_derivative_info_get="real, add more details later, #Components="!&
          !&//TRIM(NUMBER_TO_VSTRING(NUMBER_OF_COMPONENTS,"*",ERR,ERROR))
        END SELECT
      ENDDO ! dev_idx
    ENDIF !NUMBER_DERIVATIVES==1.AND.GROUP_DERIVATIVES(1)==NO_PART_DERIV

    exits("FIELD_IO_LABEL_DERIVATIVE_INFO_GET")
    RETURN
999 errorsexits("FIELD_IO_LABEL_DERIVATIVE_INFO_GET",err,error)
  END FUNCTION field_io_label_derivative_info_get

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

  FUNCTION field_io_get_field_info_label(FIELD, ERR, ERROR)
    !Argument variables
    TYPE(field_type), POINTER :: FIELD
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

    TYPE(varying_string) :: FIELD_IO_GET_FIELD_INFO_LABEL

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

    IF(.NOT.ASSOCIATED(field)) THEN
       CALL flagerror("field pointer in the input data is invalid",err,error,*999)
       GOTO 999
    ENDIF

    SELECT CASE(field%TYPE)
      CASE(field_geometric_type) !FIELD_GEOMETRIC_TYPE
        field_io_get_field_info_label="field geometric type"
      CASE(field_fibre_type)
        field_io_get_field_info_label="field fibres type"
      CASE(field_general_type)
        field_io_get_field_info_label="field general type"
      CASE(field_material_type)
        field_io_get_field_info_label="field material type"
       CASE(field_geometric_general_type)
        field_io_get_field_info_label="field geometric general type"
      CASE DEFAULT
        field_io_get_field_info_label="unknown field type"
    END SELECT

    exits("FIELD_IO_GET_FIELD_INFO_LABEL")
    RETURN
999 errorsexits("FIELD_IO_GET_FIELD_INFO_LABEL",err,error)
  END FUNCTION field_io_get_field_info_label
  !
  !================================================================================================================================
  !

  FUNCTION field_io_get_variable_info_label(COMPONENT, ERR, ERROR)
    !Argument variables
    TYPE(field_variable_component_type), POINTER :: COMPONENT
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
    TYPE(field_type), POINTER :: FIELD
    TYPE(field_variable_type), POINTER :: VARIABLE
    TYPE(varying_string) :: FIELD_IO_GET_VARIABLE_INFO_LABEL

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

    IF(.NOT.ASSOCIATED(component)) THEN
      CALL flagerror("component pointer in the input data is invalid",err,error,*999)
      GOTO 999
    ENDIF

    field=>component%FIELD_VARIABLE%FIELD
    variable=>component%FIELD_VARIABLE

    SELECT CASE(field%TYPE)
    CASE(field_geometric_type) !FIELD_GEOMETRIC_TYPE
      SELECT CASE(variable%VARIABLE_TYPE)
      CASE(field_u_variable_type)
        !coordinate system
        NULLIFY(coordinate_system)
        CALL field_coordinate_system_get(field,coordinate_system,err,error,*999)
        SELECT CASE(coordinate_system%TYPE)
          CASE(coordinate_rectangular_cartesian_type)
            field_io_get_variable_info_label="coordinates,  coordinate, rectangular cartesian"
          !CASE(COORDINATE_CYCLINDRICAL_POLAR_TYPE)
          !CASE(COORDINATE_SPHERICAL_POLAR_TYPE)
          !CASE(COORDINATE_PROLATE_SPHEROIDAL_TYPE)
          !CASE(COORDINATE_OBLATE_SPHEROIDAL_TYPE)
          CASE DEFAULT
            field_io_get_variable_info_label="unknown" !coordinates, coordinate, rectangular cartesian,
        END SELECT
      CASE(field_deludeln_variable_type)
        field_io_get_variable_info_label="Normal_derivative,  field,  normal derivative of variable"
      CASE(field_deludelt_variable_type)
        field_io_get_variable_info_label="first_time_derivative,  field,  first time derivative of variable"
      CASE(field_del2udelt2_variable_type)
        field_io_get_variable_info_label="second_time_derivative,  field,  second time derivative of variable"
      CASE DEFAULT
        field_io_get_variable_info_label="unknown_geometry,  field,  real"
      END SELECT !CASE(VARIABLE%VARIABLE_TYPE)
    CASE(field_fibre_type)
      SELECT CASE(variable%VARIABLE_TYPE)
      CASE(field_u_variable_type)
!kmith - 17.10.08: Fixing fibre field label
        !FIELD_IO_GET_VARIABLE_INFO_LABEL="fiber,  standand variable type"
        field_io_get_variable_info_label="fibres, anatomical, fibre"
!kmith - 17.10.08:
      CASE(field_deludeln_variable_type)
        field_io_get_variable_info_label="norm_der_fiber,  normal derivative of variable"
      CASE(field_deludelt_variable_type)
        field_io_get_variable_info_label="first_time_fiber,  first time derivative of variable"
      CASE(field_del2udelt2_variable_type)
        field_io_get_variable_info_label="second_time_fiber,  second time derivative of variable"
      CASE DEFAULT
        field_io_get_variable_info_label="unknown_fiber,  real"
      END SELECT !CASE(VARIABLE%VARIABLE_TYPE)
    CASE(field_general_type)
      SELECT CASE(variable%VARIABLE_TYPE)
      CASE(field_u_variable_type)
!kmith - 17.10.08: Fixing general field label
        !FIELD_IO_GET_VARIABLE_INFO_LABEL="general_variabe,  field,  string"
        field_io_get_variable_info_label="general,  field,  rectangular cartesian"
!kmith - 17.10.08:
      CASE(field_deludeln_variable_type)
        field_io_get_variable_info_label="norm_dev_variable,  field,  string"
      CASE(field_deludelt_variable_type)
        field_io_get_variable_info_label="first_time_variable,  field,  first time derivative of variable"
      CASE(field_del2udelt2_variable_type)
        field_io_get_variable_info_label="second_time_variable,  field,  second time derivative of variable"
      CASE DEFAULT
        field_io_get_variable_info_label="unknown_general,  field,  real"
      END SELECT !CASE(VARIABLE%VARIABLE_TYPE)
    CASE(field_material_type)
      SELECT CASE(variable%VARIABLE_TYPE)
      CASE(field_u_variable_type)
!kmith - 17.10.08: Fixing material field label
        !FIELD_IO_GET_VARIABLE_INFO_LABEL="material,  field,  standand variable type"
        field_io_get_variable_info_label="material,  field,  rectangular cartesian"
!kmith - 17.10.08:
      CASE(field_deludeln_variable_type)
        field_io_get_variable_info_label="normal_material,  field,  normal derivative of variable"
      CASE(field_deludelt_variable_type)
        field_io_get_variable_info_label="fist_time_material,  field,  first time derivative of variable"
      CASE(field_del2udelt2_variable_type)
        field_io_get_variable_info_label="second_time_material,  field,  second time derivative of variable"
      CASE DEFAULT
        field_io_get_variable_info_label="unknown material,  field,  real"
      END SELECT !CASE(VARIABLE%VARIABLE_TYPE)
    CASE(field_geometric_general_type)
      SELECT CASE(variable%VARIABLE_TYPE)
      CASE(field_u_variable_type)
!kmith - 17.10.08: Fixing general field label
        !FIELD_IO_GET_VARIABLE_INFO_LABEL="general_variabe,  field,  string"
        field_io_get_variable_info_label="geometric general,  field,  rectangular cartesian"
!kmith - 17.10.08:
      CASE(field_deludeln_variable_type)
        field_io_get_variable_info_label="norm_dev_variable,  field,  string"
      CASE(field_deludelt_variable_type)
        field_io_get_variable_info_label="first_time_variable,  field,  first time derivative of variable"
      CASE(field_del2udelt2_variable_type)
        field_io_get_variable_info_label="second_time_variable,  field,  second time derivative of variable"
      CASE DEFAULT
        field_io_get_variable_info_label="unknown_general,  field,  real"
      END SELECT !CASE(VARIABLE%VARIABLE_TYPE)
    CASE DEFAULT
      SELECT CASE(variable%VARIABLE_TYPE)
      CASE(field_u_variable_type)
        field_io_get_variable_info_label="unknown,  field,  unknown standand variable type"
      CASE(field_deludeln_variable_type)
        field_io_get_variable_info_label="unknown,  field,  unknown normal derivative of variable"
      CASE(field_deludelt_variable_type)
        field_io_get_variable_info_label="unknown,  field,  unknown first time derivative of variable"
      CASE(field_del2udelt2_variable_type)
        field_io_get_variable_info_label="unknown, field,  unknown second time derivative of variable"
      CASE DEFAULT
        field_io_get_variable_info_label="unknown,  field,  real"
      END SELECT !CASE(VARIABLE%VARIABLE_TYPE)
    END SELECT

    exits("FIELD_IO_GET_VARIABLE_INFO_LABEL")
    RETURN
999 errorsexits("FIELD_IO_GET_VARIABLE_INFO_LABEL",err,error)
  END FUNCTION field_io_get_variable_info_label
  !
  !================================================================================================================================
  !

  FUNCTION field_io_get_component_info_label(COMPONENT, ERR, ERROR)
    !Argument variables
    TYPE(field_variable_component_type), POINTER :: COMPONENT
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
    TYPE(field_type), POINTER :: FIELD
    TYPE(field_variable_type), POINTER :: VARIABLE
    TYPE(varying_string) :: FIELD_IO_GET_COMPONENT_INFO_LABEL

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

    IF(.NOT.ASSOCIATED(component)) THEN
       CALL flagerror("component pointer in the input data is invalid",err,error,*999)
       GOTO 999
    ENDIF

    field=>component%FIELD_VARIABLE%FIELD
    variable=>component%FIELD_VARIABLE

    SELECT CASE(field%TYPE)
    CASE(field_geometric_type) !FIELD_GEOMETRIC_TYPE
      SELECT CASE(variable%VARIABLE_TYPE)
      CASE(field_u_variable_type)
        !coordinate system
        NULLIFY(coordinate_system)
        CALL field_coordinate_system_get(field,coordinate_system,err,error,*999)
        SELECT CASE(coordinate_system%TYPE)
        CASE(coordinate_rectangular_cartesian_type)
          IF(component%COMPONENT_NUMBER==1) THEN
            field_io_get_component_info_label="x"
          ELSE IF(component%COMPONENT_NUMBER==2) THEN
            field_io_get_component_info_label="y"
          ELSE IF(component%COMPONENT_NUMBER==3) THEN
            field_io_get_component_info_label="z"
          ENDIF
        !CASE(COORDINATE_CYCLINDRICAL_POLAR_TYPE)
        !CASE(COORDINATE_SPHERICAL_POLAR_TYPE)
        !CASE(COORDINATE_PROLATE_SPHEROIDAL_TYPE)
        !CASE(COORDINATE_OBLATE_SPHEROIDAL_TYPE)
        CASE DEFAULT
          field_io_get_component_info_label=trim(number_to_vstring(component%COMPONENT_NUMBER,"*",err,error))
        END SELECT
      CASE DEFAULT
        field_io_get_component_info_label=trim(number_to_vstring(component%COMPONENT_NUMBER,"*",err,error))
      END SELECT !CASE(VARIABLE%VARIABLE_TYPE)
      CASE DEFAULT
        field_io_get_component_info_label=trim(number_to_vstring(component%COMPONENT_NUMBER,"*",err,error))
    END SELECT

    exits("FIELD_IO_GET_COMPONENT_INFO_LABEL")
    RETURN
999 errorsexits("FIELD_IO_GET_COMPONENT_INFO_LABEL",err,error)
  END FUNCTION field_io_get_component_info_label

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

  !!>Write the header of a group nodes using FORTRAIN
  !SUBROUTINE FIELD_IO_IMPORT_NODAL_GROUP_HEADER_FORTRAN(NODAL_INFO_SET, LOCAL_NODAL_NUMBER, MAX_NUM_OF_NODAL_DERIVATIVES, &
  !&my_computational_node_number, FILE_ID, ERR,ERROR, *)
  !  !Argument variables
  !  TYPE(FIELD_IO_INFO_SET), INTENT(INOUT) :: NODAL_INFO_SET  !<NODAL_INFO_SET
  !  INTEGER(INTG), INTENT(IN) :: LOCAL_NODAL_NUMBER !<LOCAL_NUMBER IN THE NODAL IO LIST
  !  INTEGER(INTG), INTENT(INOUT) :: MAX_NUM_OF_NODAL_DERIVATIVES !<MAX_NUM_OF_NODAL_DERIVATIVES
  !  INTEGER(INTG), INTENT(IN) :: my_computational_node_number !<local process number
  !  INTEGER(INTG), INTENT(IN) :: FILE_ID !< FILE ID
  !  INTEGER(INTG), INTENT(OUT) :: ERR !<The error code
  !  TYPE(VARYING_STRING), INTENT(OUT) :: ERROR !<The error string
  !  !Local Variables
  !  TYPE(FIELD_TYPE), POINTER :: field_ptr
  !  TYPE(FIELD_VARIABLE_TYPE), POINTER :: variable_ptr
  !  TYPE(DOMAIN_MAPPING_TYPE), POINTER :: DOMAIN_MAPPING_NODES !The domain mapping to calculate nodal mappings
  !  TYPE(DOMAIN_NODES_TYPE), POINTER :: DOMAIN_NODES ! domain nodes
  !  TYPE(VARYING_STRING) :: LINE, LABEL
  !  INTEGER(INTG) :: NUM_OF_FIELDS, NUM_OF_VARIABLES, NUM_OF_NODAL_DEV
  !  INTEGER(INTG) :: local_number, global_number
  !  INTEGER(INTG), POINTER :: GROUP_FIELDS(:), GROUP_VARIABLES(:), GROUP_DERIVATIVES(:)
  !  INTEGER(INTG) :: field_idx, comp_idx, comp_idx1, value_idx, var_idx, global_var_idx !dev_idx,
  !
  !  ENTERS("FIELD_IO_IMPORT_NODAL_GROUP_HEADER_FORTRAN",ERR,ERROR,*999)
  !
  !  !colllect nodal header information for IO first
  !
  !  !!get the number of this computational node from mpi pool
  !  !my_computational_node_number=COMPUTATIONAL_NODE_NUMBER_GET(ERR,ERROR)
  !  !IF(ERR/=0) GOTO 999
  !
  !  !attach the temporary pointer
  !  !tmp_components=>NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%COMPONENTS
  !
  !  !collect maximum number of nodal derivatives, number of fields and variables
  !  NUM_OF_FIELDS=0
  !  NUM_OF_VARIABLES=0
  !  MAX_NUM_OF_NODAL_DERIVATIVES=0
  !  global_number=NODAL_INFO_SET%LIST_OF_GLOBAL_NUMBER(LOCAL_NODAL_NUMBER)
  !  NULLIFY(field_ptr)
  !  NULLIFY(variable_ptr)
  !  DO comp_idx=1,NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%NUMBER_OF_COMPONENTS
  !     !calculate the number of fields
  !     IF (.NOT.ASSOCIATED(field_ptr, target=NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)% &
  !          &COMPONENTS(comp_idx)%PTR%FIELD)) THEN
  !        NUM_OF_FIELDS=NUM_OF_FIELDS+1
  !        field_ptr=>NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%COMPONENTS (comp_idx)%PTR%FIELD
  !     ENDIF
  !
  !     !calculate the number of variables
  !     IF (.NOT.ASSOCIATED(variable_ptr, target=NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)% &
  !          &COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE)) THEN
  !        NUM_OF_VARIABLES=NUM_OF_VARIABLES+1
  !        variable_ptr=>NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE
  !     ENDIF
  !
  !    !finding the local numbering through the global to local mapping
  !     DOMAIN_MAPPING_NODES=>NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%&
  !         &COMPONENTS(comp_idx)%PTR%DOMAIN%MAPPINGS%NODES
  !     !get the domain index for this variable component according to my own computional node number
  !     local_number = FindMyLocalDomainNumber( DOMAIN_MAPPING_NODES%GLOBAL_TO_LOCAL_MAP(global_number), my_computational_node_number )
  !     !use local domain information find the out the maximum number of derivatives
  !     DOMAIN_NODES=>NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%COMPONENTS(comp_idx)%PTR%DOMAIN%TOPOLOGY%NODES
  !     MAX_NUM_OF_NODAL_DERIVATIVES=MAX(DOMAIN_NODES%NODES(local_number)%NUMBER_OF_DERIVATIVES,MAX_NUM_OF_NODAL_DERIVATIVES)
  !  ENDDO !comp_idx
  !  !Allocate the memory for group of field variables
  !  ALLOCATE(GROUP_FIELDS(NUM_OF_FIELDS),STAT=ERR)
  !  IF(ERR/=0) CALL FlagError("Could not allocate temporary field buffer in IO",ERR,ERROR,*999)
  !  !Allocate the memory for group of field components
  !  ALLOCATE(GROUP_VARIABLES(NUM_OF_VARIABLES),STAT=ERR)
  !  IF(ERR/=0) CALL FlagError("Could not allocate temporary variable buffer in IO",ERR,ERROR,*999)
  !  !Allocate the memory for group of maximum number of derivatives
  !  ALLOCATE(GROUP_DERIVATIVES(MAX_NUM_OF_NODAL_DERIVATIVES),STAT=ERR)
  !  IF(ERR/=0) CALL FlagError("Could not allocate temporary derivatives buffer in IO",ERR,ERROR,*999)
  !
  !  !fill information into the group of fields and variables
  !  NUM_OF_FIELDS=0
  !  NUM_OF_VARIABLES=0
  !  NULLIFY(field_ptr)
  !  NULLIFY(variable_ptr)
  !  GROUP_FIELDS(:)=0 !the item in this arrary is the number of variables in the same field
  !  GROUP_VARIABLES(:)=0 !the item in this arrary is the number of components in the same variable
  !  DO comp_idx=1,NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%NUMBER_OF_COMPONENTS
  !     !grouping field variables and components together
  !     IF((.NOT.ASSOCIATED(field_ptr,TARGET=NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)% &
  !          &COMPONENTS(comp_idx)%PTR%FIELD)).AND.(.NOT.ASSOCIATED(variable_ptr,TARGET=NODAL_INFO_SET% &
  !          &NODAL_INFO_SET(LOCAL_NODAL_NUMBER)%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE))) THEN !different field and variables
  !        !add one new variable
  !        NUM_OF_FIELDS=NUM_OF_FIELDS+1
  !        GROUP_FIELDS(NUM_OF_FIELDS)=GROUP_FIELDS(NUM_OF_FIELDS)+1
  !        !add one new component
  !        NUM_OF_VARIABLES=NUM_OF_VARIABLES+1
  !        GROUP_VARIABLES(NUM_OF_VARIABLES)=GROUP_VARIABLES(NUM_OF_VARIABLES)+1
  !        field_ptr=>NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%COMPONENTS(comp_idx)%PTR%FIELD
  !        variable_ptr=>NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE
  !     ELSE IF (ASSOCIATED(field_ptr,TARGET=NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)% &
  !        &COMPONENTS(comp_idx)%PTR%FIELD).AND.(.NOT.ASSOCIATED(variable_ptr,TARGET=NODAL_INFO_SET%&
  !        &NODAL_INFO_SET(LOCAL_NODAL_NUMBER)%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE))) THEN !the same field and  different variables
  !        !add one new variable
  !        GROUP_FIELDS(NUM_OF_FIELDS)=GROUP_FIELDS(NUM_OF_FIELDS)+1
  !        !add one new component
  !        NUM_OF_VARIABLES=NUM_OF_VARIABLES+1
  !        GROUP_VARIABLES(NUM_OF_VARIABLES)=GROUP_VARIABLES(NUM_OF_VARIABLES)+1
  !        variable_ptr=>NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE
  !     ELSE  !different components of the same variable
  !        !add one new component
  !        GROUP_VARIABLES(NUM_OF_VARIABLES)=GROUP_VARIABLES(NUM_OF_VARIABLES)+1
  !     ENDIF !field_ptr/=NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%COMPONENTS%COMPONENTS(comp_idx)%PTR%FIELD
  !  ENDDO  !comp_idx
  !
  !  !write out the nodal header
  !  var_idx=1
  !  comp_idx=1
  !  field_idx=1
  !  value_idx=1
  !  comp_idx1=1
  !  global_var_idx=0
  !
  !  CALL FIELD_IO_FORTRAN_FILE_READ_STRING(FILE_ID, LINE, FILE_END, ERR,ERROR, *999)
  !  IF(LINE/=" "//"#Fields="//TRIM(NUMBER_TO_VSTRING(SUM(GROUP_FIELDS(1:NUM_OF_FIELDS)),"*",ERR,ERROR))) &
  !     & CALL FlagError("Fields number in the Header part do not match",ERR,ERROR,*999)
  !
  !  DO field_idx=1, NUM_OF_FIELDS
  !     !write out the field information
  !     !LABEL=FIELD_IO_GET_FIELD_INFO_LABEL(NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%COMPONENTS(comp_idx1)%PTR, FIELD_IO_FIELD_LABEL,ERR,ERROR)
  !     !IF(ERR/=0) THEN
  !     !   CALL FlagError("can not get field label",ERR,ERROR,*999)
  !     !ENDIF
  !     !CALL FIELD_IO_FORTRAN_FILE_READ_STRING(FILE_ID, LINE, FILE_END, ERR,ERROR, *999)
  !     !IF(LINE/=TRIM(NUMBER_TO_VSTRING(field_idx,"*",ERR,ERROR))//") "//TRIM(LABEL)&
  !     !&//" , #variables="//TRIM(NUMBER_TO_VSTRING(GROUP_FIELDS(field_idx),"*",ERR,ERROR)) &
  !     !CALL FlagError("Variable number in the Header part do not match",ERR,ERROR,*999)
  !
  !     DO var_idx=1, GROUP_FIELDS(field_idx)
  !        global_var_idx=global_var_idx+1
  !        !write out the field information
  !        LABEL="  "//TRIM(NUMBER_TO_VSTRING(global_var_idx,"*",ERR,ERROR))//") "&
  !        &//FIELD_IO_GET_FIELD_INFO_LABEL(NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%&
  !        &COMPONENTS(comp_idx1)%PTR, FIELD_IO_VARIABLE_LABEL,ERR,ERROR)
  !        IF(ERR/=0) THEN
  !           CALL FlagError("can not get variable label",ERR,ERROR,*999)
  !           GOTO 999
  !        ENDIF
  !        CALL FIELD_IO_FORTRAN_FILE_READ_STRING(FILE_ID, LINE, FILE_END, ERR,ERROR, *999)
  !        IF(LINE/=TRIM(LABEL)//", #Components="//TRIM(NUMBER_TO_VSTRING(GROUP_VARIABLES(global_var_idx),"*",ERR,ERROR))) &
  !           & CALL FlagError("Components number in the Header part do not match",ERR,ERROR,*999)
  !
  !
  !        DO comp_idx=1, GROUP_VARIABLES(global_var_idx)
  !           !write out the component information
  !           LABEL="   "//FIELD_IO_GET_FIELD_INFO_LABEL(NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%&
  !           &COMPONENTS(comp_idx1)%PTR, FIELD_IO_COMPONENT_LABEL,ERR,ERROR)
  !           IF(ERR/=0) THEN
  !              CALL FlagError("can not get component label",ERR,ERROR,*999)
  !              GOTO 999
  !           ENDIF
  !           LINE=TRIM(LABEL)//"."
  !
  !           !finding the local numbering through the global to local mapping
  !           DOMAIN_MAPPING_NODES=>NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%COMPONENTS(comp_idx)%PTR%&
  !              &DOMAIN%MAPPINGS%NODES
  !           !get the domain index for this variable component according to my own computional node number
  !           local_number = FindMyLocalDomainNumber( DOMAIN_MAPPING_NODES%GLOBAL_TO_LOCAL_MAP(global_number), my_computational_node_number )
  !           !use local domain information find the out the maximum number of derivatives
  !           DOMAIN_NODES=>NODAL_INFO_SET%COMPONENT_INFO_SET(LOCAL_NODAL_NUMBER)%COMPONENTS(comp_idx)%PTR%DOMAIN%TOPOLOGY%NODES
  !           !get the nodal partial derivatives
  !           NUM_OF_NODAL_DEV=DOMAIN_NODES%NODES(local_number)%NUMBER_OF_DERIVATIVES
  !           GROUP_DERIVATIVES(1:NUM_OF_NODAL_DEV)=DOMAIN_NODES%NODES(local_number)%PARTIAL_DERIVATIVE_INDEX(:)
  !           !sort  the partial derivatives
  !           CALL LIST_SORT(GROUP_DERIVATIVES(1:NUM_OF_NODAL_DEV),ERR,ERROR,*999)
  !           !get the derivative name
  !           LABEL=FIELD_IO_LABEL_DERIVATIVE_INFO_GET(GROUP_DERIVATIVES(1:NUM_OF_NODAL_DEV), NUM_OF_NODAL_DEV, &
  !              &FIELD_IO_DERIVATIVE_LABEL,ERR,ERROR)
  !           IF(ERR/=0) THEN
  !              CALL FlagError("can not get derivative label",ERR,ERROR,*999)
  !              GOTO 999
  !           ENDIF
  !           !write out the header
  !           CALL FIELD_IO_FORTRAN_FILE_READ_STRING(FILE_ID, LINE, FILE_END, ERR,ERROR, *999)
  !           !assemble the header
  !           IF(LINE/=LINE//"  Value index= "//TRIM(NUMBER_TO_VSTRING(value_idx,"*",ERR,ERROR))&
  !              &//", #Derivatives= "//TRIM(NUMBER_TO_VSTRING(NUM_OF_NODAL_DEV-1,"*",ERR,ERROR))//TRIM(LABEL)) &
  !              & CALL FlagError("Value index in the Header part do not match",ERR,ERROR,*999)
  !           !increase the component index
  !           comp_idx1=comp_idx1+1
  !           !increase the value index
  !           value_idx=value_idx+NUM_OF_NODAL_DEV
  !        ENDDO !comp_idx
  !     ENDDO !var_idx
  !  ENDDO !field_idx
  !
  !  !release temporary memory
  !  IF(ASSOCIATED(GROUP_FIELDS)) DEALLOCATE(GROUP_FIELDS)
  !  IF(ASSOCIATED(GROUP_VARIABLES)) DEALLOCATE(GROUP_VARIABLES)
  !  IF(ASSOCIATED(GROUP_DERIVATIVES)) DEALLOCATE(GROUP_DERIVATIVES)
  !
  !  EXITS("FIELD_IO_IMPORT_NODAL_GROUP_HEADER_FORTRAN")
  !  RETURN
!999 ERRORSEXITS("FIELD_IO_IMPORT_NODAL_GROUP_HEADER_FORTRAN",ERR,ERROR)
  !  RETURN 1
  !END SUBROUTINE FIELD_IO_IMPORT_NODAL_GROUP_HEADER_FORTRAN

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

  SUBROUTINE field_io_export_nodal_group_header_fortran(fieldInfoSet, global_number, MAX_NUM_OF_NODAL_DERIVATIVES, &
  &my_computational_node_number, sessionhandle, paddinginfo, err,error, *)
    !Argument variables
    TYPE(field_io_component_info_set), INTENT(IN) :: fieldInfoSet
    INTEGER(INTG), INTENT(IN) :: global_number
    INTEGER(INTG), INTENT(INOUT) :: MAX_NUM_OF_NODAL_DERIVATIVES
    INTEGER(INTG), INTENT(IN) :: my_computational_node_number
    INTEGER(INTG), INTENT(IN) :: sessionHandle
    INTEGER(INTG), ALLOCATABLE, INTENT(INOUT) :: paddingInfo(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: i, LENGTH
    CHARACTER(LEN=MAXSTRLEN) :: fvar_name
    CHARACTER(LEN=1, KIND=C_CHAR) :: cvar_name(maxstrlen+1)
    TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
    TYPE(field_type), POINTER :: field_ptr
    TYPE(field_variable_type), POINTER :: variable_ptr
    TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES ! domain nodes
    TYPE(field_variable_component_type), POINTER :: component, fieldComponent
    INTEGER(INTG), ALLOCATABLE, TARGET :: GROUP_FIELDS(:), GROUP_VARIABLES(:), GROUP_DERIVATIVES(:)
    INTEGER(INTG) :: NUM_OF_FIELDS, NUM_OF_VARIABLES, NUM_OF_NODAL_DEV
    INTEGER(INTG) :: local_number
    INTEGER(INTG) :: field_idx, comp_idx, comp_idx1, value_idx, var_idx, global_var_idx ,derivative_idx !dev_idx,
    LOGICAL :: FOUND

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

    !colllect nodal header information for IO first

    !!get the number of this computational node from mpi pool
    !my_computational_node_number=COMPUTATIONAL_NODE_NUMBER_GET(ERR,ERROR)
    !IF(ERR/=0) GOTO 999

    !attach the temporary pointer
    !tmp_components=>fieldInfoSet%COMPONENTS

    !collect maximum number of nodal derivatives, number of fields and variables
    num_of_fields=0
    num_of_variables=0
    max_num_of_nodal_derivatives=0
    NULLIFY(field_ptr)
    NULLIFY(variable_ptr)
    DO comp_idx=1,fieldinfoset%NUMBER_OF_COMPONENTS
      !calculate the number of fields
      IF (.NOT.ASSOCIATED(field_ptr, TARGET=fieldinfoset%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE%FIELD)) THEN
        num_of_fields=num_of_fields+1
        field_ptr=>fieldinfoset%COMPONENTS (comp_idx)%PTR%FIELD_VARIABLE%FIELD
      ENDIF

      !calculate the number of variables
      IF (.NOT.ASSOCIATED(variable_ptr, TARGET=fieldinfoset%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE)) THEN
        num_of_variables=num_of_variables+1
        variable_ptr=>fieldinfoset%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE
      ENDIF

      !find the local numbering
      domain_nodes=>fieldinfoset%COMPONENTS(comp_idx)%PTR%DOMAIN%TOPOLOGY%NODES
      found=.false.
      DO local_number=1,domain_nodes%NUMBER_OF_NODES
        IF( domain_nodes%NODES(local_number)%GLOBAL_NUMBER == global_number ) THEN
         found = .true.
         EXIT
        ENDIF
      ENDDO !local_number

      IF( .NOT. found ) THEN
        !Something's gone horribly wrong
        cycle
      ENDIF

      max_num_of_nodal_derivatives=max(domain_nodes%NODES(local_number)%NUMBER_OF_DERIVATIVES,max_num_of_nodal_derivatives)
    ENDDO !comp_idx

    !Allocate the memory for group of field variables
    ALLOCATE(group_fields(num_of_fields),stat=err)
    IF(err/=0) CALL flagerror("Could not allocate temporary field buffer in IO",err,error,*999)
    !Allocate the memory for group of field components
    ALLOCATE(group_variables(num_of_variables),stat=err)
    IF(err/=0) CALL flagerror("Could not allocate temporary variable buffer in IO",err,error,*999)
    !Allocate the memory for group of maximum number of derivatives
    ALLOCATE(group_derivatives(max_num_of_nodal_derivatives),stat=err)
    IF(err/=0) CALL flagerror("Could not allocate temporary derivatives buffer in IO",err,error,*999)

    !fill information into the group of fields and variables
    num_of_fields=0
    num_of_variables=0
    NULLIFY(field_ptr)
    NULLIFY(variable_ptr)
    group_fields(:)=0 !the item in this arrary is the number of variables in the same field
    group_variables(:)=0 !the item in this arrary is the number of components in the same variable
    DO comp_idx=1,fieldinfoset%NUMBER_OF_COMPONENTS
      !grouping field variables and components together
      IF((.NOT.ASSOCIATED(field_ptr,TARGET=fieldinfoset%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE%FIELD)).AND. &
          & (.NOT.ASSOCIATED(variable_ptr,TARGET=fieldinfoset%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE))) THEN !different field and variables
        num_of_fields=num_of_fields+1
        field_ptr=>fieldinfoset%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE%FIELD
      ENDIF

      IF(.NOT.ASSOCIATED(variable_ptr,TARGET=fieldinfoset%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE)) THEN !the same field and  different variables
        !add one new variable
        group_fields(num_of_fields)=group_fields(num_of_fields)+1
        !add one new component
        num_of_variables=num_of_variables+1
        variable_ptr=>fieldinfoset%COMPONENTS(comp_idx)%PTR%FIELD_VARIABLE
      ENDIF

      group_variables(num_of_variables)=group_variables(num_of_variables)+1

    ENDDO  !comp_idx

    !write out the nodal header
    var_idx=1
    comp_idx=1
    field_idx=1
    value_idx=1
    comp_idx1=1
    global_var_idx=0

    CALL reallocate( paddinginfo, fieldinfoset%NUMBER_OF_COMPONENTS + 1, "Cannot allocate padding info", err, error, *999 )

    err = fieldexport_fieldcount( sessionhandle, sum(group_fields(1:num_of_fields) ) )
    IF(err/=0) THEN
      CALL flagerror( "File write error during field export", err, error,*999 )
    ENDIF

    DO field_idx=1, num_of_fields
      DO var_idx=1, group_fields(field_idx)
        global_var_idx=global_var_idx+1

        variable_ptr=>fieldinfoset%COMPONENTS(comp_idx1)%PTR%FIELD_VARIABLE
        !write out the field information

        fvar_name = char(variable_ptr%variable_label)
        length=len_trim(fvar_name)
        DO i=1,length
          cvar_name(i)=fvar_name(i:i)
        ENDDO !i
        cvar_name(length+1)=c_null_char

        IF( variable_ptr%FIELD%TYPE == field_geometric_type .AND. &
          & variable_ptr%VARIABLE_TYPE == field_u_variable_type ) THEN
          NULLIFY(coordinate_system)
          CALL field_coordinate_system_get(variable_ptr%FIELD,coordinate_system,err,error,*999)
          err = fieldexport_coordinatevariable( sessionhandle, cvar_name, global_var_idx, &
           & coordinate_system%TYPE, variable_ptr%NUMBER_OF_COMPONENTS )
        ELSE
          err = fieldexport_variable( sessionhandle, cvar_name, global_var_idx, variable_ptr%FIELD%TYPE,  &
           & variable_ptr%VARIABLE_TYPE, &
           & variable_ptr%NUMBER_OF_COMPONENTS )
        ENDIF
        IF( err /= 0 ) THEN
          CALL flagerror( "File write error during field export", err, error,*999 )
        ENDIF

        DO comp_idx=1, variable_ptr%NUMBER_OF_COMPONENTS
          !write out the component information

          fieldcomponent => variable_ptr%COMPONENTS(comp_idx)

          IF( comp_idx1 <= fieldinfoset%NUMBER_OF_COMPONENTS ) THEN
            !It's possible to run out of node-local components before we've examined all field components.
            component => fieldinfoset%COMPONENTS(comp_idx1)%PTR
          ENDIF

          !The field component is not present at this node. Add a dummy value.
          IF(.NOT.ASSOCIATED(component,TARGET=fieldcomponent)) THEN
            paddinginfo(comp_idx1) = paddinginfo(comp_idx1) + 1
            group_derivatives(1:1) = no_part_deriv
            IF( fieldcomponent%FIELD_VARIABLE%FIELD%TYPE == field_geometric_type .AND. &
               & fieldcomponent%FIELD_VARIABLE%VARIABLE_TYPE == field_u_variable_type ) THEN
             NULLIFY(coordinate_system)
             CALL field_coordinate_system_get(variable_ptr%FIELD,coordinate_system,err,error,*999)
             err = fieldexport_coordinatederivativeindices( sessionhandle, fieldcomponent%COMPONENT_NUMBER, &
                 & coordinate_system%TYPE, 1, c_loc(group_derivatives), value_idx )
            ELSE
              err = fieldexport_derivativeindices( sessionhandle, fieldcomponent%COMPONENT_NUMBER, &
                  & variable_ptr%FIELD%TYPE, &
                  & variable_ptr%VARIABLE_TYPE, 1, c_loc(group_derivatives), value_idx )
            ENDIF

            value_idx = value_idx + 1

            err = fieldexport_endcomponent( sessionhandle )
            cycle
          ENDIF

          !use local domain information find the out the maximum number of derivatives
          domain_nodes=>component%DOMAIN%TOPOLOGY%NODES

          found=.false.
          DO local_number=1,domain_nodes%NUMBER_OF_NODES
            IF( domain_nodes%NODES(local_number)%GLOBAL_NUMBER == global_number ) THEN
              found = .true.
              EXIT
            ENDIF
          ENDDO !local_number

          IF( .NOT. found ) THEN
            err = fieldexport_endcomponent( sessionhandle )
            cycle
          ENDIF

          !get the nodal partial derivatives
          num_of_nodal_dev=domain_nodes%NODES(local_number)%NUMBER_OF_DERIVATIVES
          DO derivative_idx=1,num_of_nodal_dev
            group_derivatives(derivative_idx)=domain_nodes%NODES(local_number)%DERIVATIVES(derivative_idx)% &
                & partial_derivative_index
          ENDDO
          !sort  the partial derivatives
          CALL list_sort(group_derivatives(1:num_of_nodal_dev),err,error,*999)

          IF( component%FIELD_VARIABLE%FIELD%TYPE == field_geometric_type .AND. &
            & component%FIELD_VARIABLE%VARIABLE_TYPE == field_u_variable_type ) THEN
            NULLIFY(coordinate_system)
            CALL field_coordinate_system_get(variable_ptr%FIELD,coordinate_system,err,error,*999)
            err = fieldexport_coordinatederivativeindices( sessionhandle, component%COMPONENT_NUMBER, &
                & coordinate_system%TYPE, num_of_nodal_dev, c_loc(group_derivatives), value_idx )
          ELSE
            err = fieldexport_derivativeindices( sessionhandle, component%COMPONENT_NUMBER, &
                & variable_ptr%FIELD%TYPE, &
                & variable_ptr%VARIABLE_TYPE,num_of_nodal_dev, c_loc(group_derivatives), value_idx )
          ENDIF

          err = fieldexport_versioninfo( sessionhandle, fieldinfoset%COMPONENT_VERSIONS(comp_idx1) )
          IF(err/=0) THEN
            CALL flagerror( "Error exporting version information.", err, error,*999 )
          ENDIF
          err = fieldexport_endcomponent( sessionhandle )

          !increase the component index
          comp_idx1=comp_idx1+1
          !increase the value index
          value_idx=value_idx+num_of_nodal_dev
        ENDDO !comp_idx
      ENDDO !var_idx
    ENDDO !field_idx

    !release temporary memory
    CALL checked_deallocate( group_fields )
    CALL checked_deallocate( group_variables )
    CALL checked_deallocate( group_derivatives )

    exits("FIELD_IO_EXPORT_NODAL_GROUP_HEADER_FORTRAN")
    RETURN
999 errorsexits("FIELD_IO_EXPORT_NODAL_GROUP_HEADER_FORTRAN",err,error)
    RETURN 1
  END SUBROUTINE field_io_export_nodal_group_header_fortran

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

  SUBROUTINE field_io_export_nodes_into_local_file(NODAL_INFO_SET, NAME, my_computational_node_number,ERR, ERROR, *)
    !the reason that my_computational_node_number is used in the argument is for future extension
    !Argument variables
    TYPE(field_io_info_set), INTENT(INOUT):: NODAL_INFO_SET
    TYPE(varying_string), INTENT(IN) :: NAME
    INTEGER(INTG), INTENT(IN):: my_computational_node_number
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: FILE_NAME !the prefix name of file.
    TYPE(field_variable_component_type), POINTER :: COMPONENT !the prefix name of file.
    TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES ! domain nodes
    INTEGER(INTG) :: local_number, global_number, sessionHandle, paddingCount,DERIVATIVE_INDEXES(part_deriv_s4_s4_s4)
    INTEGER(INTG), ALLOCATABLE :: paddingInfo(:)
    INTEGER(INTG) :: nn, comp_idx,  dev_idx, version_idx, NUM_OF_NODAL_DEV, MAX_NUM_OF_NODAL_DERIVATIVES, total_nodal_values
    INTEGER(INTG) :: NUMBER_VERSIONS, MAX_NUMBER_VERSIONS
    INTEGER(INTG), POINTER :: GEOMETRIC_PARAMETERS_INTG(:)
    LOGICAL :: FOUND
    REAL(C_DOUBLE), ALLOCATABLE, TARGET :: NODAL_BUFFER(:), TOTAL_NODAL_BUFFER(:)
    REAL(DP), POINTER :: GEOMETRIC_PARAMETERS_DP(:)
    !INTEGER(INTG), POINTER :: GEOMETRIC_PARAMETERS_INTG(:)
    REAL(DP) :: padding(1),VALUE

    padding(1) = 1.23456789


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

    !get my own computianal node number--be careful the rank of process in the MPI pool
    !is not necessarily equal to numbering of computional node, so use method COMPUTATIONAL_NODE_NUMBER_GET
    !will be a secured way to get the number
    file_name=name//".part"//trim(number_to_vstring(my_computational_node_number,"*",err,error))//".exnode"
    max_num_of_nodal_derivatives=0

    IF(.NOT.ALLOCATED(nodal_info_set%COMPONENT_INFO_SET)) THEN
      CALL flagerror("the nodal information set in input is invalid",err,error,*999)
    ENDIF

    IF(.NOT.ALLOCATED(nodal_info_set%LIST_OF_GLOBAL_NUMBER)) THEN
      CALL flagerror("the nodal global information set is not associated with any numbering list",err,error,*999)
    ENDIF

    IF(nodal_info_set%NUMBER_OF_ENTRIES==0) THEN
      CALL flagerror("the nodal information set does not contain any nodes",err,error,*999)
    ENDIF

    IF(nodal_info_set%COMPONENT_INFO_SET(1)%PTR%SAME_HEADER) THEN
      CALL flagerror("the first header flag of nodal information set should be false",err,error,*999)
    ENDIF

    err = fieldexport_opensession( export_type_file, char(file_name)//c_null_char, sessionhandle )
    IF(err/=0) THEN
      CALL flagerror( "Cannot open file export session", err, error,*999 )
    ENDIF

    IF(ASSOCIATED(nodal_info_set%FIELDS%REGION)) THEN
      err = fieldexport_group( sessionhandle, char(nodal_info_set%FIELDS%REGION%LABEL)//c_null_char )
    ELSE
      IF(ASSOCIATED(nodal_info_set%FIELDS%INTERFACE)) THEN
       err = fieldexport_group( sessionhandle, char(nodal_info_set%FIELDS%INTERFACE%LABEL)//c_null_char )
     ELSE
        CALL flagerror("Fields region or interface is not associated.",err,error,*999)
      ENDIF
    ENDIF
    IF(err/=0) THEN
      CALL flagerror( "Cannot write group name to nodes file", err, error,*999 )
    ENDIF

    DO nn=1, nodal_info_set%NUMBER_OF_ENTRIES
      global_number=nodal_info_set%LIST_OF_GLOBAL_NUMBER(nn)

      IF(.NOT.nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%SAME_HEADER) THEN
        !write out the nodal header

        CALL field_io_export_nodal_group_header_fortran(nodal_info_set%COMPONENT_INFO_SET(nn)%PTR, &
          & global_number, max_num_of_nodal_derivatives, my_computational_node_number, sessionhandle, &
          & paddinginfo, err,error,*999)
        max_number_versions = maxval(nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENT_VERSIONS)
        !value_idx=value_idx-1 !the len of NODAL_BUFFER
        !checking: whether need to allocate temporary memory for Io writing
        IF(ALLOCATED(nodal_buffer)) THEN
          IF(SIZE(nodal_buffer)<max_num_of_nodal_derivatives*max_number_versions) THEN
            CALL reallocate( nodal_buffer, max_num_of_nodal_derivatives*max_number_versions, &
              & "Could not allocate temporary nodal buffer in IO writing", err, error, *999 )
          ENDIF
        ELSE
          CALL reallocate( nodal_buffer, max_num_of_nodal_derivatives*max_number_versions, &
            & "Could not allocate temporary nodal buffer in IO writing", err, error, *999 )
        ENDIF
      ENDIF !NODAL_INFO_SET%COMPONENT_INFO_SET(nn)%SAME_HEADER==.FALSE.

      !write out the components' values of this node in this domain
      total_nodal_values = 0
      CALL checked_deallocate( total_nodal_buffer )
      DO comp_idx=1,nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS
        component => nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS(comp_idx)%PTR
        number_versions = nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENT_VERSIONS(comp_idx)
        domain_nodes=>component%DOMAIN%TOPOLOGY%NODES
        found=.false.
        DO local_number=1,domain_nodes%NUMBER_OF_NODES
          IF( domain_nodes%NODES(local_number)%GLOBAL_NUMBER == global_number ) THEN
            found = .true.
            EXIT
          ENDIF
        ENDDO !local_number

        IF(.NOT. found) THEN
          cycle
        ENDIF

        DO paddingcount = 1, paddinginfo( comp_idx )
          num_of_nodal_dev = 1
          nodal_buffer(1) = padding(1)

          CALL grow_array( total_nodal_buffer, num_of_nodal_dev, "Insufficient memory during I/O", err, error, *999 )
          total_nodal_buffer(total_nodal_values+1:total_nodal_values+num_of_nodal_dev) = nodal_buffer(1:num_of_nodal_dev)
          total_nodal_values = total_nodal_values + num_of_nodal_dev

          err = fieldexport_nodevalues( sessionhandle, domain_nodes%NODES(local_number)%USER_NUMBER, num_of_nodal_dev, &
            & c_loc(nodal_buffer) )
          IF(err/=0) THEN
            CALL flagerror( "Cannot write group name to nodes file", err, error,*999 )
          ENDIF
        ENDDO !paddingCount


        SELECT CASE(component%FIELD_VARIABLE%DATA_TYPE)
        CASE(field_intg_type)
          NULLIFY(geometric_parameters_intg)
          CALL field_parameter_set_data_get(component%FIELD_VARIABLE%FIELD,component%FIELD_VARIABLE%VARIABLE_TYPE, &
            & field_values_set_type,geometric_parameters_intg,err,error,*999)
        CASE(field_dp_type)
          NULLIFY(geometric_parameters_dp)
          CALL field_parameter_set_data_get(component%FIELD_VARIABLE%FIELD,component%FIELD_VARIABLE%VARIABLE_TYPE, &
            & field_values_set_type,geometric_parameters_dp,err,error,*999)
        CASE DEFAULT
          CALL flagerror("Not implemented.",err,error,*999)
        END SELECT
        !get the nodal partial derivatives
        num_of_nodal_dev=domain_nodes%NODES(local_number)%NUMBER_OF_DERIVATIVES

        !Record the dof-index of each derivative (if it is present)
        derivative_indexes = -1
        DO dev_idx=1, num_of_nodal_dev
          derivative_indexes( domain_nodes%NODES(local_number)%DERIVATIVES(dev_idx)%PARTIAL_DERIVATIVE_INDEX ) = dev_idx
        ENDDO

        !Output the dofs, sorted according to derivative index
        !Loop over versions outside of derivatives, as cmgui treats versions differently
        num_of_nodal_dev = 0
        DO version_idx=1, number_versions
          DO dev_idx=1, SIZE(derivative_indexes)
            IF( derivative_indexes( dev_idx ) == -1 ) THEN
              cycle
            ENDIF

            num_of_nodal_dev = num_of_nodal_dev + 1

            SELECT CASE(component%FIELD_VARIABLE%DATA_TYPE)
            CASE(field_intg_type)
              IF(component%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(local_number)% &
                  & derivatives(derivative_indexes(dev_idx))%NUMBER_OF_VERSIONS < version_idx) THEN
                !If the number of versions for this derivative isn't equal to the maximum number of versions for the
                !component, then fill the rest of the version data with the first version
                VALUE=REAL(GEOMETRIC_PARAMETERS_INTG( COMPONENT%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(local_number)% &
                  & DERIVATIVES(DERIVATIVE_INDEXES(dev_idx))%VERSIONS(1) ) ,DP)
              ELSE
                VALUE=REAL(GEOMETRIC_PARAMETERS_INTG( COMPONENT%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(local_number)% &
                  & DERIVATIVES(DERIVATIVE_INDEXES(dev_idx))%VERSIONS(version_idx) ) ,DP)
              ENDIF
            CASE(field_dp_type)
              !Default to version 1 of each node derivative
              IF(component%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(local_number)% &
                  & derivatives(derivative_indexes(dev_idx))%NUMBER_OF_VERSIONS < version_idx) THEN
                VALUE=geometric_parameters_dp( component%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(local_number)% &
                  & derivatives(derivative_indexes(dev_idx))%VERSIONS(1) )
              ELSE
                VALUE=geometric_parameters_dp( component%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(local_number)% &
                  & derivatives(derivative_indexes(dev_idx))%VERSIONS(version_idx) )
              ENDIF
            CASE DEFAULT
              CALL flagerror("Not implemented.",err,error,*999)
            END SELECT
            nodal_buffer( num_of_nodal_dev ) = VALUE
          ENDDO !dev_idx
        ENDDO !verion_idx

        CALL grow_array( total_nodal_buffer, num_of_nodal_dev, "Insufficient memory during I/O", err, error, *999 )
        total_nodal_buffer(total_nodal_values+1:total_nodal_values+num_of_nodal_dev) = nodal_buffer(1:num_of_nodal_dev)
        total_nodal_values = total_nodal_values + num_of_nodal_dev

        !TEMPORARY
        err = fieldexport_nodevalues( sessionhandle, domain_nodes%NODES(local_number)%USER_NUMBER, num_of_nodal_dev, &
          & c_loc(nodal_buffer) )
        IF(err/=0) THEN
          CALL flagerror( "Cannot write group name to nodes file", err, error,*999 )
        ENDIF
      ENDDO !comp_idx

      !Note that paddingInfo's size is one more than the component count
      DO paddingcount = 1, paddinginfo( nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS + 1 )
        num_of_nodal_dev = 1
        nodal_buffer(1) = padding(1)

        CALL grow_array( total_nodal_buffer, num_of_nodal_dev, "Insufficient memory during I/O", err, error, *999 )
        total_nodal_buffer(total_nodal_values+1:total_nodal_values+num_of_nodal_dev) = nodal_buffer(1:num_of_nodal_dev)
        total_nodal_values = total_nodal_values + num_of_nodal_dev

        err = fieldexport_nodevalues( sessionhandle, domain_nodes%NODES(local_number)%USER_NUMBER, num_of_nodal_dev, &
          & c_loc(nodal_buffer) )
        IF(err/=0) THEN
          CALL flagerror( "Cannot write group name to nodes file", err, error,*999 )
        ENDIF
      ENDDO !paddingCount

      !REINSTATE FOR HDF5 BLOCK DATA WRITES
      !       ERR = FieldExport_NodeValues( sessionHandle, DOMAIN_NODES%NODES(local_number)%USER_NUMBER, &
      !         & total_nodal_values, C_LOC(TOTAL_NODAL_BUFFER) )
      !       IF(ERR/=0) THEN
      !         CALL FlagError( "Cannot write group name to nodes file", ERR, ERROR,*999 )
      !       ENDIF


    ENDDO !nn

    err = fieldexport_closesession( sessionhandle )
    IF(err/=0) THEN
      CALL flagerror( "Cannot write group name to nodes file", err, error,*999 )
    ENDIF

    !release the temporary memory
    CALL checked_deallocate( nodal_buffer )
    CALL checked_deallocate( total_nodal_buffer )
    IF(ASSOCIATED(geometric_parameters_dp)) NULLIFY(geometric_parameters_dp)
    IF(ASSOCIATED(geometric_parameters_intg)) NULLIFY(geometric_parameters_intg)

    exits("FIELD_IO_EXPORT_NODES_INTO_LOCAL_FILE")
    RETURN
999 errorsexits("FIELD_IO_EXPORT_NODES_INTO_LOCAL_FILE",err,error)
    RETURN 1
  END SUBROUTINE field_io_export_nodes_into_local_file

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

  SUBROUTINE field_io_fortran_file_read_string(FILE_ID, STRING_DATA, FILE_END, ERR,ERROR,*)

    !Argument variables
    TYPE(varying_string), INTENT(INOUT) :: STRING_DATA
    INTEGER(INTG), INTENT(IN) :: FILE_ID
    LOGICAL, INTENT(INOUT) :: FILE_END
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    CHARACTER (LEN=MAXSTRLEN) :: TEMP_STR
    INTEGER(INTG) :: LEN_OF_DATA, IOS

    string_data=remove(string_data, 1, len(string_data))

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

    READ(file_id, "(A)", iostat=ios) temp_str

    IF(ios>=0)  THEN
      file_end=.false.
    ELSE
      file_end=.true.
    ENDIF

    string_data=trim(temp_str)
    len_of_data=len(string_data)

    IF(len_of_data==0) THEN
      CALL flagerror("leng of string is zero",err,error,*999)
    ENDIF

    exits("FIELD_IO_FORTRAN_FILE_READ_STRING")
    RETURN
999 errorsexits("FIELD_IO_FORTRAN_FILE_READ_STRING",err,error)

    exits("FIELD_IO_FORTRAN_FILE_READ_STRING")
    RETURN 1
  END SUBROUTINE field_io_fortran_file_read_string

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

  SUBROUTINE field_io_fortran_file_read_dp(FILE_ID, REAL_DATA, LEN_OF_DATA, FILE_END, ERR,ERROR,*)

    !Argument variables
    REAL(DP), INTENT(OUT) :: REAL_DATA(:)
    INTEGER(INTG), INTENT(IN) :: FILE_ID
    INTEGER(INTG), INTENT(IN) :: LEN_OF_DATA
    LOGICAL, INTENT(INOUT) :: FILE_END ! file end
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: DP_FMT
    INTEGER(INTG) :: IOS

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

    dp_fmt="("//trim(number_to_vstring(len_of_data,"*",err,error))//"ES)"
    READ(file_id, char(dp_fmt), iostat=ios) real_data(1:len_of_data)

    IF(ios>=0)  THEN
      file_end=.false.
    ELSE
      file_end=.true.
    ENDIF

    exits("FIELD_IO_FORTRAN_FILE_READ_DP")
    RETURN
999 errorsexits("FIELD_IO_FORTRAN_FILE_READ_DP",err,error)
    RETURN 1
  END SUBROUTINE field_io_fortran_file_read_dp

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

  SUBROUTINE field_io_fortran_file_write_dp(FILE_ID, REAL_DATA, LEN_OF_DATA, ERR,ERROR,*)

    !Argument variables
    REAL(DP), INTENT(IN) :: REAL_DATA(:)
    INTEGER(INTG), INTENT(IN) :: FILE_ID
    INTEGER(INTG), INTENT(IN) :: LEN_OF_DATA
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    !DP_FMT="(ES"//TRIM(NUMBER_TO_VSTRING(LEN_OF_DATA,"*",ERR,ERROR))//".0)"
    !WRITE(FILE_ID, CHAR(DP_FMT)) REAL_DATA(1:LEN_OF_DATA)
    WRITE(file_id,*) real_data(1:len_of_data)

    exits("FIELD_IO_FORTRAN_FILE_WRITE_DP")
    RETURN
999 errorsexits("FIELD_IO_FORTRAN_FILE_WRITE_DP",err,error)
    RETURN 1
  END SUBROUTINE field_io_fortran_file_write_dp

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

  SUBROUTINE field_io_fortran_file_read_intg(FILE_ID, INTG_DATA, LEN_OF_DATA, ERR,ERROR,*)

    !Argument variables
    INTEGER(INTG), INTENT(OUT) :: INTG_DATA(:)
    INTEGER(INTG), INTENT(IN) :: FILE_ID
    INTEGER(INTG), INTENT(IN) :: LEN_OF_DATA
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: DP_FMT

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

    dp_fmt="("//trim(number_to_vstring(len_of_data,"*",err,error))//"I)"
    READ(file_id, char(dp_fmt)) intg_data(1:len_of_data)

    exits("FIELD_IO_FORTRAN_FILE_READ_INTG")
    RETURN
999 errorsexits("FIELD_IO_FORTRAN_FILE_READ_INTG",err,error)
    RETURN 1
  END SUBROUTINE field_io_fortran_file_read_intg

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

  SUBROUTINE field_io_fortran_file_write_intg(FILE_ID, INTG_DATA, LEN_OF_DATA, ERR,ERROR,*)

    !Argument variables
    INTEGER(INTG), INTENT(IN) :: INTG_DATA(:)
    INTEGER(INTG), INTENT(IN) :: FILE_ID
    INTEGER(INTG), INTENT(IN) :: LEN_OF_DATA
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: DP_FMT

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

    dp_fmt="(I"//trim(number_to_vstring(len_of_data,"*",err,error))//")"
    WRITE(file_id, char(dp_fmt)) intg_data(1:len_of_data)

    exits("FIELD_IO_FORTRAN_FILE_WRITE_INTG")
    RETURN
999 errorsexits("FIELD_IO_FORTRAN_FILE_WRITE_INTG",err,error)
    RETURN 1
  END SUBROUTINE field_io_fortran_file_write_intg

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

  SUBROUTINE field_io_fortran_file_open(FILE_ID, FILE_NAME, FILE_STATUS, ERR,ERROR,*)

    !Argument variables
    TYPE(varying_string), INTENT(INOUT) :: FILE_NAME
    TYPE(varying_string), INTENT(IN) :: FILE_STATUS
    INTEGER(INTG), INTENT(INOUT) :: FILE_ID
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    !CALL WRITE_STRING(DIAGNOSTIC_OUTPUT_TYPE,"OPEN FILE",ERR,ERROR,*999)

    OPEN(unit=file_id, file=char(file_name), status=char(file_status), form="FORMATTED", err=999)


    exits("FIELD_IO_FORTRAN_FILE_OPEN")
    RETURN
999 errorsexits("FIELD_IO_FORTRAN_FILE_OPEN",err,error)
    RETURN 1
  END SUBROUTINE field_io_fortran_file_open

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

  SUBROUTINE field_io_fortran_file_close(FILE_ID, ERR,ERROR, *)
    !Argument variables
    INTEGER(INTG), INTENT(INOUT) :: FILE_ID
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    !CALL WRITE_STRING(DIAGNOSTIC_OUTPUT_TYPE,"CLOSE FILE",ERR,ERROR,*999)

    CLOSE(unit=file_id, err=999)

    exits("FIELD_IO_FORTRAN_FILE_CLOSE")
    RETURN
999 errorsexits("FIELD_IO_FORTRAN_FILE_CLOSE",err,error)
    RETURN 1
  END SUBROUTINE field_io_fortran_file_close

  SUBROUTINE string_to_muti_integers_vs(STRING, NUMBER_OF_INTEGERS, INTG_DATA, ERR, ERROR, *)

    !#### Function: STRING_TO_INTEGER_VS
    !###  Type: INTEGER(INTG)
    !###  Description:
    !###    Converts a varying string representation of a number to an integer.
    !###  Parent-function: STRING_TO_INTEGER

    !Argument variables
    TYPE(varying_string), INTENT(IN) :: STRING
    INTEGER(INTG), INTENT(INOUT) :: INTG_DATA(:)
    INTEGER(INTG), INTENT(IN) :: NUMBER_OF_INTEGERS
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_STRING, LOCAL_STRING1
    INTEGER(INTG) :: idx, pos

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

!!TODO: remove dependance on LOCAL_STRING

    local_string=string

    DO idx=1,number_of_integers-1
      local_string=adjustl(local_string)
      local_string=trim(local_string)
      pos=index(local_string, " ")
      local_string1=extract(local_string, 1, pos-1)
      intg_data(idx)=string_to_integer(local_string1, err, error)
      local_string=remove(local_string,1,pos)
    ENDDO
    local_string=adjustl(local_string)
    local_string=trim(local_string)
    intg_data(idx)=string_to_integer(local_string, err, error)

    exits("STRING_TO_MUTI_INTEGERS_VS")
    RETURN
999 errorsexits("STRING_TO_MUTI_INTEGERS_VS",err,error)
    RETURN 1
  END SUBROUTINE string_to_muti_integers_vs

  SUBROUTINE string_to_muti_reals_vs(STRING, NUMBER_OF_REALS, REAL_DATA, POSITION, ERR, ERROR, *)

    !#### Function: STRING_TO_INTEGER_VS
    !###  Type: INTEGER(INTG)
    !###  Description:
    !###    Converts a varying string representation of a number to an integer.
    !###  Parent-function: STRING_TO_INTEGER

    !Argument variables
    TYPE(varying_string), INTENT(IN) :: STRING
    REAL(DP), INTENT(INOUT) :: REAL_DATA(:)
    INTEGER(INTG), INTENT(IN) :: NUMBER_OF_REALS
    INTEGER(INTG), INTENT(IN) :: POSITION
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_STRING, LOCAL_STRING1
    !CHARACTER(256) :: CHAR_BUFF
    INTEGER(INTG) :: idx, pos

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

!!TODO: remove dependance on LOCAL_STRING

    local_string=string

    DO idx=1,number_of_reals-1
      local_string=adjustl(local_string)
      local_string=trim(local_string)
      pos=index(local_string, " ")
      local_string1=extract(local_string, 1, pos-1)
      !CHAR_BUFF=CHAR(LOCAL_STRING1)
      !READ(CHAR_BUFF,"(ES)",IOSTAT=ERR,ERR=999) REAL_DATA(idx+POSITION-1)
      real_data(idx+position-1)=string_to_double(local_string, err, error)
      local_string=remove(local_string,1,pos)
    ENDDO
    local_string=adjustl(local_string)
    local_string=trim(local_string)
    real_data(idx+position-1)=string_to_double(local_string, err, error)
    !CHAR_BUFF=CHAR(LOCAL_STRING)
    !READ(CHAR_BUFF,"(ES)",IOSTAT=ERR,ERR=999) REAL_DATA(idx+POSITION-1)
    !REAL_DATA(idx+POSITION-1)=STRING_TO_DOUBLE(LOCAL_STRING, ERR, ERROR)

    exits("STRING_TO_MUTI_REALS_VS")
    RETURN
999 errorsexits("STRING_TO_MUTI_REALS_VS",err,error)
    RETURN 1
  END SUBROUTINE string_to_muti_reals_vs


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

  SUBROUTINE fieldio_nodelinfosetattachlocalprocess(NODAL_INFO_SET, FIELDS, my_computational_node_number, ERR,ERROR,*)
    !Argument variables
    TYPE(field_io_info_set), INTENT(INOUT):: NODAL_INFO_SET
    TYPE(fields_type), POINTER ::FIELDS
    INTEGER(INTG), INTENT(IN):: my_computational_node_number
    INTEGER(INTG), INTENT(OUT):: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    !Local Variables
    LOGICAL :: ININTERFACE,INREGION
    TYPE(field_type), POINTER :: FIELD
    TYPE(domain_nodes_type), POINTER:: DOMAIN_NODES !nodes in local domain
    TYPE(field_variable_type), POINTER:: FIELD_VARIABLE !field variable
    INTEGER(INTG) :: field_idx, var_idx, component_idx, deriv_idx, np, nn, num_field !temporary variable
    INTEGER(INTG) :: MAX_NUMBER_VERSIONS
    LOGICAL :: foundNewNode
    TYPE(varying_string) :: LOCAL_ERROR

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

    !validate the input data
    inregion=.false.
    ininterface=.false.
    IF(ASSOCIATED(fields%REGION)) THEN
      inregion=.true.
    ELSE
      IF(ASSOCIATED(fields%INTERFACE)) THEN
        ininterface=.true.
      ELSE
        CALL flagerror("Fields is not associated with a region or interface.",err,error,*999)
      ENDIF
    ENDIF

    IF(inregion) THEN
      !checking whether the list of fields in the same region
      DO num_field =1, fields%NUMBER_OF_FIELDS
        IF(.NOT.ASSOCIATED(fields%FIELDS(num_field)%PTR)) THEN
          local_error ="No. "//trim(number_to_vstring(num_field,"*",err,error))// &
            & " field handle in fields list is invalid"
          CALL flagerror(local_error,err,error,*999)
        ENDIF

        IF( num_field == 1 ) THEN
          cycle
        ENDIF

        IF(fields%FIELDS(num_field-1)%PTR%REGION%USER_NUMBER/=fields%FIELDS(num_field)%PTR%REGION%USER_NUMBER) THEN
          local_error = "No. "//trim(number_to_vstring(num_field-1,"*",err,error))//" and "// &
            & trim(number_to_vstring(num_field,"*",err,error))//" fields are not in the same region"
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ENDDO
    ELSE
      !checking whether the list of fields in the same interface
      DO num_field =1, fields%NUMBER_OF_FIELDS
        IF(.NOT.ASSOCIATED(fields%FIELDS(num_field)%PTR)) THEN
          local_error ="No. "//trim(number_to_vstring(num_field,"*",err,error))// &
            & " field handle in fields list is invalid"
          CALL flagerror(local_error,err,error,*999)
        ENDIF

        IF( num_field == 1 ) THEN
          cycle
        ENDIF

        IF(fields%FIELDS(num_field-1)%PTR%INTERFACE%USER_NUMBER/= &
          & fields%FIELDS(num_field)%PTR%INTERFACE%USER_NUMBER) THEN
          local_error = "No. "//trim(number_to_vstring(num_field-1,"*",err,error))//" and "// &
            & trim(number_to_vstring(num_field,"*",err,error))//" fields are not in the same interface."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ENDDO
    ENDIF

    !checking whether the list of fields are using the same decomposition
    !IF(.NOT.ASSOCIATED(DECOMPOSITION))
    !  CALL FlagError("decomposition method is not vakid",ERR,ERROR,*999)
    !ENDIF
    !DO num_field =1, FIELDS%NUMBER_OF_FIELDS
    !  IF(FIELDS%FIELDS(num_field)%PTR%DECOMPOSITION/=DECOMPOSITION)
    !    LOCAL_ERROR ="No. "//TRIM(NUMBER_TO_VSTRING(num_field,"*",ERR,ERROR)) //" field "&
    !    & //" uses different decomposition method with the specified decomposition method,"//&
    !     & "which is not supported currently, ask Heye for more details"
    !    CALL FlagError(LOCAL_ERROR,ERR,ERROR,*999)
    !  ENDIF
    !ENDDO

    nodal_info_set%FIELDS=>fields

    !attach local process to local nodal information set. In current opencmiss system,
    !each local process owns it local nodal information, so all we need to do is to fill the nodal
    !information set with nodal information of local process
    IF( ( nodal_info_set%NUMBER_OF_ENTRIES > 0 ) .OR. &
      & ALLOCATED( nodal_info_set%COMPONENT_INFO_SET ) ) THEN
      CALL flagerror("nodal information set is not initialized properly, call start method first",err,error,*999)
    ENDIF

    DO field_idx = 1, nodal_info_set%FIELDS%NUMBER_OF_FIELDS
      field => nodal_info_set%FIELDS%FIELDS(field_idx)%PTR
      IF( .NOT.ALLOCATED( field%VARIABLES ) ) THEN
        cycle
      ENDIF

      DO var_idx = 1, field%NUMBER_OF_VARIABLES
        field_variable => field%VARIABLES( var_idx )

        DO component_idx = 1, field_variable%NUMBER_OF_COMPONENTS
          IF( .NOT.ASSOCIATED( field_variable%COMPONENTS( component_idx )%DOMAIN%TOPOLOGY%NODES ) ) THEN
            cycle
          ENDIF

          domain_nodes => field_variable%COMPONENTS( component_idx )%DOMAIN%TOPOLOGY%NODES
          !TODO This is an order n-squared algorithm. It would be faster to build a list of
          !non-unique nodes, sort it, then remove duplicates.

          !Add this domain's node indexes to LIST_OF_GLOBAL_NUMBER if it's not present.
          DO np = 1, domain_nodes%NUMBER_OF_NODES
            foundnewnode = .true.
            DO nn = 1,nodal_info_set%NUMBER_OF_ENTRIES
              IF( nodal_info_set%LIST_OF_GLOBAL_NUMBER( nn ) == domain_nodes%NODES( np )%GLOBAL_NUMBER ) THEN
                foundnewnode = .false.
                EXIT
              ENDIF
            ENDDO !nn=1

            IF( foundnewnode ) THEN
              !TODO This code re-allocates new memory once per node. Reduce memory-thrashing by adding memory in chunks.
              CALL grow_array( nodal_info_set%LIST_OF_GLOBAL_NUMBER, 1, "Could not allocate buffer in IO", err, error, *999 )
              nodal_info_set%LIST_OF_GLOBAL_NUMBER(nodal_info_set%NUMBER_OF_ENTRIES+1)= domain_nodes%NODES(np)%GLOBAL_NUMBER
              nodal_info_set%NUMBER_OF_ENTRIES=nodal_info_set%NUMBER_OF_ENTRIES+1
            ENDIF
          ENDDO !np
        ENDDO !component_idx
      ENDDO !var_idx
    ENDDO !field_idx

    !allocate the nodal information set and initialize them
    ALLOCATE( nodal_info_set%COMPONENT_INFO_SET( nodal_info_set%NUMBER_OF_ENTRIES ), stat = err )
    IF( err /= 0 ) CALL flagerror( "Could not allocate nodal information set", err, error, *999)

    DO nn = 1, nodal_info_set%NUMBER_OF_ENTRIES
      ALLOCATE( nodal_info_set%COMPONENT_INFO_SET(nn)%PTR )
      nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%SAME_HEADER = .false.
      nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS = 0
      CALL checked_deallocate( nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS )
      CALL checked_deallocate( nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENT_VERSIONS )
    ENDDO

    !collect nodal information from local process
    DO field_idx = 1, nodal_info_set%FIELDS%NUMBER_OF_FIELDS
      field => nodal_info_set%FIELDS%FIELDS( field_idx )%PTR
      IF( .NOT.ALLOCATED(field%VARIABLES) ) THEN
        cycle
      ENDIF

      DO var_idx=1, field%NUMBER_OF_VARIABLES
        field_variable => field%VARIABLES( var_idx )
        DO component_idx = 1, field_variable%NUMBER_OF_COMPONENTS
          IF( field_variable%COMPONENTS( component_idx )%INTERPOLATION_TYPE /= field_node_based_interpolation ) THEN
            cycle
          ENDIF

          IF( .NOT.ASSOCIATED( field_variable%COMPONENTS( component_idx )%DOMAIN%TOPOLOGY%NODES ) ) THEN
            cycle
          ENDIF

          domain_nodes => field_variable%COMPONENTS( component_idx )%DOMAIN%TOPOLOGY%NODES

          DO np = 1, domain_nodes%NUMBER_OF_NODES
            DO nn = 1, nodal_info_set%NUMBER_OF_ENTRIES
              IF( nodal_info_set%LIST_OF_GLOBAL_NUMBER( nn ) == domain_nodes%NODES( np )%GLOBAL_NUMBER ) THEN
                ! Old CMISS and cmgui treat versions differently, the version loop is outside the derivative
                ! loop rather than the other way around as in OpenCMISS, so we have to have the same number of
                ! versions for all derivatives
                max_number_versions = 1
                DO deriv_idx=1,domain_nodes%NODES( np )%NUMBER_OF_DERIVATIVES
                  max_number_versions = max(max_number_versions, &
                    & domain_nodes%NODES(np)%DERIVATIVES(deriv_idx)%numberOfVersions)
                END DO
                !allocate variable component memory
                CALL grow_array( nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS, 1, &
                  & "Could not allocate temporary buffer in IO", err, error, *999 )
                nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS(nodal_info_set%COMPONENT_INFO_SET(nn)%PTR &
                  & %NUMBER_OF_COMPONENTS+1)%PTR=>field%VARIABLES( var_idx )%COMPONENTS( component_idx )
                nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS = &
                  & nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS+1
                CALL grow_array( nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENT_VERSIONS, 1, &
                  & "Could not allocate temporary buffer in IO", err, error, *999 )
                nodal_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENT_VERSIONS( &
                  & nodal_info_set%COMPONENT_INFO_SET(nn)%PTR &
                  & %NUMBER_OF_COMPONENTS) = max_number_versions
                EXIT
              ENDIF
            ENDDO
          ENDDO !np
        ENDDO !component_idx
      ENDDO !var_idx
    ENDDO !field_idx

    exits("FieldIO_NodelInfoSetAttachLocalProcess")
    RETURN
999 errorsexits("FieldIO_NodelInfoSetAttachLocalProcess",err,error)
    RETURN 1
  END SUBROUTINE fieldio_nodelinfosetattachlocalprocess

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

  SUBROUTINE field_io_info_set_initialise( LOCAL_PROCESS_INFO_SET, ERR, ERROR, * )
    !Argument variables
    TYPE(field_io_info_set), INTENT(INOUT) :: LOCAL_PROCESS_INFO_SET
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: nn, ncomp  !temporary variable

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

    IF(ASSOCIATED(local_process_info_set%FIELDS)) THEN
       NULLIFY(local_process_info_set%FIELDS)
    ENDIF
    IF(ALLOCATED(local_process_info_set%COMPONENT_INFO_SET)) THEN
      DO nn=1, local_process_info_set%NUMBER_OF_ENTRIES
        IF(ALLOCATED(local_process_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS)) THEN
          DO ncomp=1, local_process_info_set%COMPONENT_INFO_SET(nn)%PTR%NUMBER_OF_COMPONENTS
            NULLIFY(local_process_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS(ncomp)%PTR)
          ENDDO
          CALL checked_deallocate( local_process_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENTS )
          CALL checked_deallocate( local_process_info_set%COMPONENT_INFO_SET(nn)%PTR%COMPONENT_VERSIONS )
          DEALLOCATE( local_process_info_set%COMPONENT_INFO_SET(nn)%PTR )
        ENDIF
      ENDDO
      DEALLOCATE(local_process_info_set%COMPONENT_INFO_SET)
    ENDIF

    local_process_info_set%NUMBER_OF_ENTRIES=0
    CALL checked_deallocate( local_process_info_set%LIST_OF_GLOBAL_NUMBER )

    exits("FIELD_IO_INFO_SET_INITIALISE")
    RETURN
999 errorsexits("FIELD_IO_INFO_SET_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE field_io_info_set_initialise

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

  SUBROUTINE field_io_nodes_export(FIELDS, FILE_NAME, METHOD, ERR,ERROR,*)
    !Argument variables
    TYPE(fields_type), POINTER :: FIELDS
    TYPE(varying_string), INTENT(IN) :: FILE_NAME
    TYPE(varying_string), INTENT(IN):: METHOD
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(field_io_info_set) :: NODAL_INFO_SET
    INTEGER(INTG):: my_computational_node_number
    INTEGER(INTG):: computational_node_numbers

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

    !Get the number of computational nodes
    computational_node_numbers=computational_nodes_number_get(err,error)
    IF(err/=0) GOTO 999
    !Get my computational node number
    my_computational_node_number=computational_node_number_get(err,error)
    IF(err/=0) GOTO 999
    IF(method=="FORTRAN") THEN
      CALL field_io_info_set_initialise(nodal_info_set, err,error,*999)
      CALL fieldio_nodelinfosetattachlocalprocess(nodal_info_set, fields, my_computational_node_number, err,error,*999)
      CALL field_io_nodal_info_set_sort(nodal_info_set, my_computational_node_number, err,error,*999)
      CALL field_io_export_nodes_into_local_file(nodal_info_set, file_name, my_computational_node_number, &
          & err, error, *999)
      CALL field_io_info_set_initialise(nodal_info_set, err,error,*999)
    ELSE IF(method=="MPIIO") THEN
      CALL flagerror("MPI IO has not been implemented yet!",err,error,*999)
    ELSE
      CALL flagerror("Unknown method!",err,error,*999)
    ENDIF

    exits("FIELD_IO_NODES_EXPORT")
    RETURN
999 errorsexits("FIELD_IO_NODES_EXPORT",err,error)
    RETURN 1
  END SUBROUTINE field_io_nodes_export

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

  SUBROUTINE field_io_elements_export(FIELDS, FILE_NAME, METHOD,ERR,ERROR,*)
  !checking the input data for IO and initialize the nodal information set
  !the following items will be checked: the region (the same?), all the pointer(valid?)
  !in this version, different decomposition method will be allowed for the list of field variables(but still in the same region)
  !even the each process has exactly the same nodal information and each process will write out exactly the same data
  !because CMGui can read the same data for several times

    !Argument variables
    TYPE(fields_type), POINTER :: FIELDS
    TYPE(varying_string), INTENT(IN) :: FILE_NAME
    TYPE(varying_string), INTENT(IN):: METHOD
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(field_io_info_set) :: LOCAL_PROCESS_ELEMENTAL_INFO_SET
    INTEGER(INTG):: my_computational_node_number
    INTEGER(INTG):: computational_node_numbers

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

    !Get the number of computational nodes
    computational_node_numbers=computational_nodes_number_get(err,error)
    IF(err/=0) GOTO 999
    !Get my computational node number
    my_computational_node_number=computational_node_number_get(err,error)
    IF(err/=0) GOTO 999
    IF(method=="FORTRAN") THEN
      CALL field_io_info_set_initialise( local_process_elemental_info_set, err, error, *999 )
      CALL fieldio_elementalinfosetattachlocalprocess( local_process_elemental_info_set, fields, err, error, *999 )
      CALL field_io_elemental_info_set_sort(local_process_elemental_info_set, my_computational_node_number, err,error,*999)
      CALL field_io_export_elements_into_local_file(local_process_elemental_info_set, file_name, my_computational_node_number, &
          & err, error, *999)
      CALL field_io_info_set_initialise(local_process_elemental_info_set, err,error,*999)
    ELSE IF(method=="MPIIO") THEN
      CALL flagerror("MPI IO has not been implemented yet",err,error,*999)
    ELSE
      CALL flagerror("Unknown method!",err,error,*999)
    ENDIF

    exits("FIELD_IO_ELEMENTS_EXPORT")
    RETURN
999 errorsexits("FIELD_IO_ELEMENTS_EXPORT",err,error)
    RETURN 1
  END SUBROUTINE field_io_elements_export

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

END MODULE field_io_routines