anchor.diff:    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)
  +  ! Read field headers
  +  E_IO = VTK_VAR_XML_HEADER_READ(var_location,varname,NC_NN,NCOMP,fmt,type,data,offs,np,rf)
  +
  +  if(E_IO == 0) then
  +    allocate(var(NC_NN*NCOMP), stat=E_IO)
  +
  +    ! Read field data
  +    select case(vtk(rf)%f)
  +      case(ascii)
           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
  +
  +      case(binary)
  +! Not implemented
  +
  +      case(raw)
  +        if(E_IO == 0) then
  +          call get_int(buffer=s_buffer, attrib='offset', val=offs, 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
   
  -    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
  +    endselect
  +  endif
     !---------------------------------------------------------------------------------------------------------------------------------
     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).
  +  !< Function for reading field of scalar variable (I1P, 1D array).
     !---------------------------------------------------------------------------------------------------------------------------------
     implicit none
     character(*),              intent(IN)  :: var_location ! location of variables: CELL for cell-centered, NODE for node-centered
  
positive.code:   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
  
negative.code:           trim(adjustlt(Upper_Case(type)))/='INT32') then
          E_IO = -1_I4P
        else
          read(data, fmt=*, iostat=E_IO) offset 
          ! get appended array cell_type
          allocate(cell_type(NC), stat=E_IO)
          E_IO = search(from=pos,inside='Cells', to_find='DataArray', with_attribute='Name', of_value='Types', &
  