anchor.diff:              call get_char(buffer=s_buffer, attrib='WholeExtent', val=aux, E_IO=E_IO)
               if(E_IO == 0) then
                 read(aux,*) rn
  -              if(present(nx1)) nx1 = rn(1)
  -              if(present(nx2)) nx2 = rn(2)
  -              if(present(ny1)) ny1 = rn(3)
  -              if(present(ny2)) ny2 = rn(4)
  -              if(present(nz1)) nz1 = rn(5)
  -              if(present(nz2)) nz2 = rn(6)
  +              if(present(nx1)) nx1 = rn(1); if(present(nx2)) nx2 = rn(2)
  +              if(present(ny1)) ny1 = rn(3); if(present(ny2)) ny2 = rn(4)
  +              if(present(nz1)) nz1 = rn(5); if(present(nz2)) nz2 = rn(6)
               endif
           end select
   
  
positive.code:         allocate(cell_type(NC), stat=E_IO)
          E_IO = search(from=pos,inside='Cells', to_find='DataArray', with_attribute='Name', of_value='Types', &
                        buffer=s_buffer, content=data)
          call get_int(buffer=s_buffer,attrib='offset', val=offs, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='format', val=fmt, E_IO=E_IO)
          if (trim(adjustlt(Upper_Case(fmt)))/='ASCII') then
            E_IO = -1_I4P
          else
            read(data, fmt=*, iostat=E_IO) cell_type !get appended array connect
  
negative.code:   !---------------------------------------------------------------------------------------------------------------------------------
  end function
  
  !  module procedure VTK_VAR_XML_SCAL_1DA_R8,VTK_VAR_XML_SCAL_3DA_R8, & ! real(R8P)    scalar    1D/3D array
  !                   VTK_VAR_XML_SCAL_1DA_R4,VTK_VAR_XML_SCAL_3DA_R4, & ! real(R4P)    scalar    1D/3D array
  !                   VTK_VAR_XML_SCAL_1DA_I8,VTK_VAR_XML_SCAL_3DA_I8, & ! integer(I8P) scalar    1D/3D array
  !                   VTK_VAR_XML_SCAL_1DA_I4,VTK_VAR_XML_SCAL_3DA_I4, & ! integer(I4P) scalar    1D/3D array
  !                   VTK_VAR_XML_SCAL_1DA_I2,VTK_VAR_XML_SCAL_3DA_I2, & ! integer(I2P) scalar    1D/3D array
  !                   VTK_VAR_XML_SCAL_1DA_I1,VTK_VAR_XML_SCAL_3DA_I1, & ! integer(I1P) scalar    1D/3D array
  !                   VTK_VAR_XML_VECT_1DA_R8,VTK_VAR_XML_VECT_3DA_R8, & ! real(R8P)    vectorial 1D/3D arrays
  !                   VTK_VAR_XML_VECT_1DA_R4,VTK_VAR_XML_VECT_3DA_R4, & ! real(R4P)    vectorial 1D/3D arrays
  !                   VTK_VAR_XML_VECT_1DA_I8,VTK_VAR_XML_VECT_3DA_I8, & ! integer(I8P) vectorial 1D/3D arrays
  !                   VTK_VAR_XML_VECT_1DA_I4,VTK_VAR_XML_VECT_3DA_I4, & ! integer(I4P) vectorial 1D/3D arrays
  !                   VTK_VAR_XML_VECT_1DA_I2,VTK_VAR_XML_VECT_3DA_I2, & ! integer(I2P) vectorial 1D/3D arrays
  !                   VTK_VAR_XML_VECT_1DA_I1,VTK_VAR_XML_VECT_3DA_I1, & ! integer(I1P) vectorial 1D/3D arrays
  !                   VTK_VAR_XML_LIST_1DA_R8,VTK_VAR_XML_LIST_3DA_R8, & ! real(R8P)    list      1D/3D array
  !                   VTK_VAR_XML_LIST_1DA_R4,VTK_VAR_XML_LIST_3DA_R4, & ! real(R4P)    list      1D/3D array
  !                   VTK_VAR_XML_LIST_1DA_I8,VTK_VAR_XML_LIST_3DA_I8, & ! integer(I4P) list      1D/3D array
  !                   VTK_VAR_XML_LIST_1DA_I4,VTK_VAR_XML_LIST_3DA_I4, & ! integer(I4P) list      1D/3D array
  !                   VTK_VAR_XML_LIST_1DA_I2,VTK_VAR_XML_LIST_3DA_I2, & ! integer(I2P) list      1D/3D array
  !                   VTK_VAR_XML_LIST_1DA_I1,VTK_VAR_XML_LIST_3DA_I1    ! integer(I1P) list      1D/3D array
  !-----------------------------------------------------------------------------------------------------------------------------------
    function VTK_VAR_XML_SCAL_1DA_R8_READ(var_location,varname,NC_NN,NCOMP,var,npiece,cf) result(E_IO)
    !---------------------------------------------------------------------------------------------------------------------------------
    !< Function for reading field of scalar variable (R8P, 1D array).
    !---------------------------------------------------------------------------------------------------------------------------------
    implicit none
    character(*),           intent(IN)  :: var_location ! location of variables: CELL for cell-centered, NODE for node-centered
    character(*),           intent(IN)  :: varname      ! variable name
    integer(I4P),           intent(OUT) :: NC_NN        ! number of cells or nodes
    integer(I4P),           intent(OUT) :: NCOMP        ! number of components
    real(R8P), allocatable, intent(OUT) :: var(:)       ! variable to be saved
    integer(I4P), optional, intent(IN)  :: npiece       ! Number of the piece to read (by default: 1)
    integer(I4P), optional, intent(IN)  :: cf           !< Current file index (for concurrent files IO).
    integer(I4P)                        :: E_IO         ! Input/Output inquiring flag: $0$ if IO is done, $> 0$ if IO is not done
    integer(I4P)                        :: rf           !< Real file index.
    character(len=:), allocatable       :: s_buffer     !< Buffer string.
    character(len=:), allocatable       :: fmt
    character(len=:), allocatable       :: type
    character(len=:), allocatable       :: data
    integer(I4P)                        :: np, offs, pos
    integer(I4P)                        :: N_Byte
    !---------------------------------------------------------------------------------------------------------------------------------
  
    !---------------------------------------------------------------------------------------------------------------------------------
    E_IO = -1_I4P
    rf = f
    if (present(cf)) then
      rf = cf ; f = cf
    endif
    np = 1_I4P; if (present(npiece)) np = npiece
  
    select case(vtk(rf)%f)
      case(ascii)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside=trim(vtk(rf)%topology), to_find='Piece', repeat=np, buffer=s_buffer)
        inquire(unit=vtk(rf)%u, pos=pos)
        select case(trim(Upper_case(var_location)))
        case('NODE')
          call get_int(buffer=s_buffer, attrib='NumberOfPoints', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='PointData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer, content=data)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        case('CELL')
          call get_int(buffer=s_buffer, attrib='NumberOfCells', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='CellData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer, content=data)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        end select
        if(E_IO == 0) then
          call get_int(buffer=s_buffer, attrib='offset', val=offs, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='format', val=fmt, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='type', val=type, E_IO=E_IO)
          if (trim(adjustlt(Upper_Case(fmt)))/='ASCII' .or. &
              trim(adjustlt(Upper_Case(type)))/='FLOAT64') then
            E_IO = -1_I4P
          else
            read(data, fmt=*, iostat=E_IO) var
          endif
        endif
  
      case(binary)
  ! Not implemented
  
      case(raw)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside=trim(vtk(rf)%topology), to_find='Piece', repeat=np, buffer=s_buffer)
        inquire(unit=vtk(rf)%u, pos=pos)
        select case(trim(Upper_case(var_location)))
        case('NODE')
          call get_int(buffer=s_buffer, attrib='NumberOfPoints', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='PointData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        case('CELL')
          call get_int(buffer=s_buffer, attrib='NumberOfCells', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='CellData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        end select
        if(E_IO == 0) then
          call get_int(buffer=s_buffer, attrib='offset', val=offs, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='format', val=fmt, E_IO=E_IO)
          if (trim(adjustlt(Upper_Case(fmt)))/='APPENDED' .or. &
              trim(adjustlt(Upper_Case(type)))/='FLOAT64') then
            E_IO = -1_I4P
          else
            read(unit=vtk(rf)%u, iostat=E_IO, pos=vtk(rf)%ioffset+offs) N_Byte, var
          endif
        endif
  
    endselect
    !---------------------------------------------------------------------------------------------------------------------------------
    end function
  
  
    function VTK_VAR_XML_SCAL_1DA_R4_READ(var_location,varname,NC_NN,NCOMP,var,npiece,cf) result(E_IO)
    !---------------------------------------------------------------------------------------------------------------------------------
    !< Function for reading field of scalar variable (R8P, 1D array).
    !---------------------------------------------------------------------------------------------------------------------------------
    implicit none
    character(*),           intent(IN)  :: var_location ! location of variables: CELL for cell-centered, NODE for node-centered
    character(*),           intent(IN)  :: varname      ! variable name
    integer(I4P),           intent(OUT) :: NC_NN        ! number of cells or nodes
    integer(I4P),           intent(OUT) :: NCOMP        ! number of components
    real(R4P), allocatable, intent(OUT) :: var(:)       ! variable to be saved
    integer(I4P), optional, intent(IN)  :: npiece       ! Number of the piece to read (by default: 1)
    integer(I4P), optional, intent(IN)  :: cf           !< Current file index (for concurrent files IO).
    integer(I4P)                        :: E_IO         ! Input/Output inquiring flag: $0$ if IO is done, $> 0$ if IO is not done
    integer(I4P)                        :: rf           !< Real file index.
    character(len=:), allocatable       :: s_buffer     !< Buffer string.
    character(len=:), allocatable       :: fmt
    character(len=:), allocatable       :: type
    character(len=:), allocatable       :: data
    integer(I4P)                        :: np, offs, pos
    integer(I4P)                        :: N_Byte
    !---------------------------------------------------------------------------------------------------------------------------------
  
    !---------------------------------------------------------------------------------------------------------------------------------
    E_IO = -1_I4P
    rf = f
    if (present(cf)) then
      rf = cf ; f = cf
    endif
    np = 1_I4P; if (present(npiece)) np = npiece
  
    select case(vtk(rf)%f)
      case(ascii)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside=trim(vtk(rf)%topology), to_find='Piece', repeat=np, buffer=s_buffer)
        inquire(unit=vtk(rf)%u, pos=pos)
        select case(trim(Upper_case(var_location)))
        case('NODE')
          call get_int(buffer=s_buffer, attrib='NumberOfPoints', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='PointData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer, content=data)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        case('CELL')
          call get_int(buffer=s_buffer, attrib='NumberOfCells', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='CellData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer, content=data)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        end select
        if(E_IO == 0) then
          call get_int(buffer=s_buffer, attrib='offset', val=offs, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='format', val=fmt, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='type', val=type, E_IO=E_IO)
          if (trim(adjustlt(Upper_Case(fmt)))/='ASCII' .or. &
              trim(adjustlt(Upper_Case(type)))/='FLOAT32') then
            E_IO = -1_I4P
          else
            read(data, fmt=*, iostat=E_IO) var
          endif
        endif
  
      case(binary)
  ! Not implemented
  
      case(raw)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside=trim(vtk(rf)%topology), to_find='Piece', repeat=np, buffer=s_buffer)
        inquire(unit=vtk(rf)%u, pos=pos)
        select case(trim(Upper_case(var_location)))
        case('NODE')
          call get_int(buffer=s_buffer, attrib='NumberOfPoints', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='PointData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        case('CELL')
          call get_int(buffer=s_buffer, attrib='NumberOfCells', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='CellData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        end select
        if(E_IO == 0) then
          call get_int(buffer=s_buffer, attrib='offset', val=offs, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='format', val=fmt, E_IO=E_IO)
          if (trim(adjustlt(Upper_Case(fmt)))/='APPENDED' .or. &
              trim(adjustlt(Upper_Case(type)))/='FLOAT32') then
            E_IO = -1_I4P
          else
            read(unit=vtk(rf)%u, iostat=E_IO, pos=vtk(rf)%ioffset+offs) N_Byte, var
          endif
        endif
  
    endselect
    !---------------------------------------------------------------------------------------------------------------------------------
    end function
  
  
    function VTK_VAR_XML_SCAL_1DA_I8_READ(var_location,varname,NC_NN,NCOMP,var,npiece,cf) result(E_IO)
    !---------------------------------------------------------------------------------------------------------------------------------
    !< Function for reading field of scalar variable (R8P, 1D array).
    !---------------------------------------------------------------------------------------------------------------------------------
    implicit none
    character(*),              intent(IN)  :: var_location ! location of variables: CELL for cell-centered, NODE for node-centered
    character(*),              intent(IN)  :: varname      ! variable name
    integer(I4P),              intent(OUT) :: NC_NN        ! number of cells or nodes
    integer(I4P),              intent(OUT) :: NCOMP        ! number of components
    integer(I8P), allocatable, intent(OUT) :: var(:)       ! variable to be saved
    integer(I4P), optional,    intent(IN)  :: npiece       ! Number of the piece to read (by default: 1)
    integer(I4P), optional,    intent(IN)  :: cf           !< Current file index (for concurrent files IO).
    integer(I4P)                           :: E_IO         ! Input/Output inquiring flag: $0$ if IO is done, $> 0$ if IO is not done
    integer(I4P)                           :: rf           !< Real file index.
    character(len=:), allocatable          :: s_buffer     !< Buffer string.
    character(len=:), allocatable          :: fmt
    character(len=:), allocatable          :: type
    character(len=:), allocatable          :: data
    integer(I4P)                           :: np, offs, pos
    integer(I4P)                           :: N_Byte
    !---------------------------------------------------------------------------------------------------------------------------------
  
    !---------------------------------------------------------------------------------------------------------------------------------
    E_IO = -1_I4P
    rf = f
    if (present(cf)) then
      rf = cf ; f = cf
    endif
    np = 1_I4P; if (present(npiece)) np = npiece
  
    select case(vtk(rf)%f)
      case(ascii)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside=trim(vtk(rf)%topology), to_find='Piece', repeat=np, buffer=s_buffer)
        inquire(unit=vtk(rf)%u, pos=pos)
        select case(trim(Upper_case(var_location)))
        case('NODE')
          call get_int(buffer=s_buffer, attrib='NumberOfPoints', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='PointData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer, content=data)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        case('CELL')
          call get_int(buffer=s_buffer, attrib='NumberOfCells', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='CellData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer, content=data)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        end select
        if(E_IO == 0) then
          call get_int(buffer=s_buffer, attrib='offset', val=offs, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='format', val=fmt, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='type', val=type, E_IO=E_IO)
          if (trim(adjustlt(Upper_Case(fmt)))/='ASCII' .or. &
              trim(adjustlt(Upper_Case(type)))/='INT64') then
            E_IO = -1_I4P
          else
            read(data, fmt=*, iostat=E_IO) var
          endif
        endif
  
      case(binary)
  ! Not implemented
  
      case(raw)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside=trim(vtk(rf)%topology), to_find='Piece', repeat=np, buffer=s_buffer)
        inquire(unit=vtk(rf)%u, pos=pos)
        select case(trim(Upper_case(var_location)))
        case('NODE')
          call get_int(buffer=s_buffer, attrib='NumberOfPoints', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='PointData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        case('CELL')
          call get_int(buffer=s_buffer, attrib='NumberOfCells', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='CellData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        end select
        if(E_IO == 0) then
          call get_int(buffer=s_buffer, attrib='offset', val=offs, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='format', val=fmt, E_IO=E_IO)
          if (trim(adjustlt(Upper_Case(fmt)))/='APPENDED' .or. &
              trim(adjustlt(Upper_Case(type)))/='INT64') then
            E_IO = -1_I4P
          else
            read(unit=vtk(rf)%u, iostat=E_IO, pos=vtk(rf)%ioffset+offs) N_Byte, var
          endif
        endif
  
    endselect
    !---------------------------------------------------------------------------------------------------------------------------------
    end function
  
    function VTK_VAR_XML_SCAL_1DA_I4_READ(var_location,varname,NC_NN,NCOMP,var,npiece,cf) result(E_IO)
    !---------------------------------------------------------------------------------------------------------------------------------
    !< Function for reading field of scalar variable (R8P, 1D array).
    !---------------------------------------------------------------------------------------------------------------------------------
    implicit none
    character(*),              intent(IN)  :: var_location ! location of variables: CELL for cell-centered, NODE for node-centered
    character(*),              intent(IN)  :: varname      ! variable name
    integer(I4P),              intent(OUT) :: NC_NN        ! number of cells or nodes
    integer(I4P),              intent(OUT) :: NCOMP        ! number of components
    integer(I4P), allocatable, intent(OUT) :: var(:)       ! variable to be saved
    integer(I4P), optional,    intent(IN)  :: npiece       ! Number of the piece to read (by default: 1)
    integer(I4P), optional,    intent(IN)  :: cf           !< Current file index (for concurrent files IO).
    integer(I4P)                           :: E_IO         ! Input/Output inquiring flag: $0$ if IO is done, $> 0$ if IO is not done
    integer(I4P)                           :: rf           !< Real file index.
    character(len=:), allocatable          :: s_buffer     !< Buffer string.
    character(len=:), allocatable          :: fmt
    character(len=:), allocatable          :: type
    character(len=:), allocatable          :: data
    integer(I4P)                           :: np, offs, pos
    integer(I4P)                           :: N_Byte
    !---------------------------------------------------------------------------------------------------------------------------------
  
    !---------------------------------------------------------------------------------------------------------------------------------
    E_IO = -1_I4P
    rf = f
    if (present(cf)) then
      rf = cf ; f = cf
    endif
    np = 1_I4P; if (present(npiece)) np = npiece
  
    select case(vtk(rf)%f)
      case(ascii)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside=trim(vtk(rf)%topology), to_find='Piece', repeat=np, buffer=s_buffer)
        inquire(unit=vtk(rf)%u, pos=pos)
        select case(trim(Upper_case(var_location)))
        case('NODE')
          call get_int(buffer=s_buffer, attrib='NumberOfPoints', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='PointData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer, content=data)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        case('CELL')
          call get_int(buffer=s_buffer, attrib='NumberOfCells', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='CellData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer, content=data)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        end select
        if(E_IO == 0) then
          call get_int(buffer=s_buffer, attrib='offset', val=offs, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='format', val=fmt, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='type', val=type, E_IO=E_IO)
          if (trim(adjustlt(Upper_Case(fmt)))/='ASCII' .or. &
              trim(adjustlt(Upper_Case(type)))/='INT32') then
            E_IO = -1_I4P
          else
            read(data, fmt=*, iostat=E_IO) var
          endif
        endif
  
      case(binary)
  ! Not implemented
  
      case(raw)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside=trim(vtk(rf)%topology), to_find='Piece', repeat=np, buffer=s_buffer)
        inquire(unit=vtk(rf)%u, pos=pos)
        select case(trim(Upper_case(var_location)))
        case('NODE')
          call get_int(buffer=s_buffer, attrib='NumberOfPoints', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='PointData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        case('CELL')
          call get_int(buffer=s_buffer, attrib='NumberOfCells', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='CellData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        end select
        if(E_IO == 0) then
          call get_int(buffer=s_buffer, attrib='offset', val=offs, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='format', val=fmt, E_IO=E_IO)
          if (trim(adjustlt(Upper_Case(fmt)))/='APPENDED' .or. &
              trim(adjustlt(Upper_Case(type)))/='INT32') then
            E_IO = -1_I4P
          else
            read(unit=vtk(rf)%u, iostat=E_IO, pos=vtk(rf)%ioffset+offs) N_Byte, var
          endif
        endif
  
    endselect
    !---------------------------------------------------------------------------------------------------------------------------------
    end function
  
  
    function VTK_VAR_XML_SCAL_1DA_I2_READ(var_location,varname,NC_NN,NCOMP,var,npiece,cf) result(E_IO)
    !---------------------------------------------------------------------------------------------------------------------------------
    !< Function for reading field of scalar variable (R8P, 1D array).
    !---------------------------------------------------------------------------------------------------------------------------------
    implicit none
    character(*),              intent(IN)  :: var_location ! location of variables: CELL for cell-centered, NODE for node-centered
    character(*),              intent(IN)  :: varname      ! variable name
    integer(I4P),              intent(OUT) :: NC_NN        ! number of cells or nodes
    integer(I4P),              intent(OUT) :: NCOMP        ! number of components
    integer(I2P), allocatable, intent(OUT) :: var(:)       ! variable to be saved
    integer(I4P), optional,    intent(IN)  :: npiece       ! Number of the piece to read (by default: 1)
    integer(I4P), optional,    intent(IN)  :: cf           !< Current file index (for concurrent files IO).
    integer(I4P)                           :: E_IO         ! Input/Output inquiring flag: $0$ if IO is done, $> 0$ if IO is not done
    integer(I4P)                           :: rf           !< Real file index.
    character(len=:), allocatable          :: s_buffer     !< Buffer string.
    character(len=:), allocatable          :: fmt
    character(len=:), allocatable          :: type
    character(len=:), allocatable          :: data
    integer(I4P)                           :: np, offs, pos
    integer(I4P)                           :: N_Byte
    !---------------------------------------------------------------------------------------------------------------------------------
  
    !---------------------------------------------------------------------------------------------------------------------------------
    E_IO = -1_I4P
    rf = f
    if (present(cf)) then
      rf = cf ; f = cf
    endif
    np = 1_I4P; if (present(npiece)) np = npiece
  
    select case(vtk(rf)%f)
      case(ascii)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside=trim(vtk(rf)%topology), to_find='Piece', repeat=np, buffer=s_buffer)
        inquire(unit=vtk(rf)%u, pos=pos)
        select case(trim(Upper_case(var_location)))
        case('NODE')
          call get_int(buffer=s_buffer, attrib='NumberOfPoints', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='PointData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer, content=data)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        case('CELL')
          call get_int(buffer=s_buffer, attrib='NumberOfCells', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='CellData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer, content=data)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        end select
        if(E_IO == 0) then
          call get_int(buffer=s_buffer, attrib='offset', val=offs, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='format', val=fmt, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='type', val=type, E_IO=E_IO)
          if (trim(adjustlt(Upper_Case(fmt)))/='ASCII' .or. &
              trim(adjustlt(Upper_Case(type)))/='INT16') then
            E_IO = -1_I4P
          else
            read(data, fmt=*, iostat=E_IO) var
          endif
        endif
  
      case(binary)
  ! Not implemented
  
      case(raw)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside=trim(vtk(rf)%topology), to_find='Piece', repeat=np, buffer=s_buffer)
        inquire(unit=vtk(rf)%u, pos=pos)
        select case(trim(Upper_case(var_location)))
        case('NODE')
          call get_int(buffer=s_buffer, attrib='NumberOfPoints', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='PointData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        case('CELL')
          call get_int(buffer=s_buffer, attrib='NumberOfCells', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='CellData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        end select
        if(E_IO == 0) then
          call get_int(buffer=s_buffer, attrib='offset', val=offs, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='format', val=fmt, E_IO=E_IO)
          if (trim(adjustlt(Upper_Case(fmt)))/='APPENDED' .or. &
              trim(adjustlt(Upper_Case(type)))/='INT16') then
            E_IO = -1_I4P
          else
            read(unit=vtk(rf)%u, iostat=E_IO, pos=vtk(rf)%ioffset+offs) N_Byte, var
          endif
        endif
  
    endselect
    !---------------------------------------------------------------------------------------------------------------------------------
    end function
  
  
    function VTK_VAR_XML_SCAL_1DA_I1_READ(var_location,varname,NC_NN,NCOMP,var,npiece,cf) result(E_IO)
    !---------------------------------------------------------------------------------------------------------------------------------
    !< Function for reading field of scalar variable (R8P, 1D array).
    !---------------------------------------------------------------------------------------------------------------------------------
    implicit none
    character(*),              intent(IN)  :: var_location ! location of variables: CELL for cell-centered, NODE for node-centered
    character(*),              intent(IN)  :: varname      ! variable name
    integer(I4P),              intent(OUT) :: NC_NN        ! number of cells or nodes
    integer(I4P),              intent(OUT) :: NCOMP        ! number of components
    integer(I1P), allocatable, intent(OUT) :: var(:)       ! variable to be saved
    integer(I4P), optional,    intent(IN)  :: npiece       ! Number of the piece to read (by default: 1)
    integer(I4P), optional,    intent(IN)  :: cf           !< Current file index (for concurrent files IO).
    integer(I4P)                           :: E_IO         ! Input/Output inquiring flag: $0$ if IO is done, $> 0$ if IO is not done
    integer(I4P)                           :: rf           !< Real file index.
    character(len=:), allocatable          :: s_buffer     !< Buffer string.
    character(len=:), allocatable          :: fmt
    character(len=:), allocatable          :: type
    character(len=:), allocatable          :: data
    integer(I4P)                           :: np, offs, pos
    integer(I4P)                           :: N_Byte
    !---------------------------------------------------------------------------------------------------------------------------------
  
    !---------------------------------------------------------------------------------------------------------------------------------
    E_IO = -1_I4P
    rf = f
    if (present(cf)) then
      rf = cf ; f = cf
    endif
    np = 1_I4P; if (present(npiece)) np = npiece
  
    select case(vtk(rf)%f)
      case(ascii)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside=trim(vtk(rf)%topology), to_find='Piece', repeat=np, buffer=s_buffer)
        inquire(unit=vtk(rf)%u, pos=pos)
        select case(trim(Upper_case(var_location)))
        case('NODE')
          call get_int(buffer=s_buffer, attrib='NumberOfPoints', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='PointData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer, content=data)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        case('CELL')
          call get_int(buffer=s_buffer, attrib='NumberOfCells', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='CellData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer, content=data)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        end select
        if(E_IO == 0) then
          call get_int(buffer=s_buffer, attrib='offset', val=offs, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='format', val=fmt, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='type', val=type, E_IO=E_IO)
          if (trim(adjustlt(Upper_Case(fmt)))/='ASCII' .or. &
              trim(adjustlt(Upper_Case(type)))/='INT8') then
            E_IO = -1_I4P
          else
            read(data, fmt=*, iostat=E_IO) var
          endif
        endif
  
      case(binary)
  ! Not implemented
  
      case(raw)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside=trim(vtk(rf)%topology), to_find='Piece', repeat=np, buffer=s_buffer)
        inquire(unit=vtk(rf)%u, pos=pos)
        select case(trim(Upper_case(var_location)))
        case('NODE')
          call get_int(buffer=s_buffer, attrib='NumberOfPoints', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='PointData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        case('CELL')
          call get_int(buffer=s_buffer, attrib='NumberOfCells', val=NC_NN, E_IO=E_IO)
          E_IO = search(from=pos, inside='CellData', to_find='DataArray', with_attribute='Name', of_value=varname, &
                        buffer=s_buffer)
          if(E_IO == 0) then
            call get_int(buffer=s_buffer, attrib='NumberOfComponents', val=NCOMP, E_IO=E_IO)
            allocate(var(NC_NN*NCOMP), stat=E_IO)
          endif
        end select
        if(E_IO == 0) then
          call get_int(buffer=s_buffer, attrib='offset', val=offs, E_IO=E_IO)
          call get_char(buffer=s_buffer, attrib='format', val=fmt, E_IO=E_IO)
          if (trim(adjustlt(Upper_Case(fmt)))/='APPENDED' .or. &
              trim(adjustlt(Upper_Case(type)))/='INT8') then
            E_IO = -1_I4P
          else
            read(unit=vtk(rf)%u, iostat=E_IO, pos=vtk(rf)%ioffset+offs) N_Byte, var
          endif
        endif
  
    endselect
    !---------------------------------------------------------------------------------------------------------------------------------
    end function
  
  
  
    function VTK_END_XML_READ(cf) result(E_IO)
    !---------------------------------------------------------------------------------------------------------------------------------
  