anchor.diff:    endif
     np = 1_I4P; if (present(npiece)) np = npiece
   
  -  ! Read field headers
  -  E_IO = VTK_VAR_XML_HEADER_READ(var_location,varname,NC_NN,NCOMP,&
  -              nx1=nx1,nx2=nx2,ny1=ny1,ny2=ny2,nz1=nz1,nz2=nz2, &
  -              fmt=fmt,type=type,data=data,offs=offs,npiece=np,cf=rf)
  -
  -  if(E_IO == 0) then
  -    allocate(var(nx1:nx2,ny1:ny2,nz1:nz2), 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
  -          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
  +  select case(trim(vtk(rf)%topology))
  +    case('RectilinearGrid','StructuredGrid')
  +      ! Read field headers
  +      E_IO = VTK_VAR_XML_HEADER_READ(var_location,varname,NC_NN,NCOMP,&
  +                  nx1=nx1,nx2=nx2,ny1=ny1,ny2=ny2,nz1=nz1,nz2=nz2, &
  +                  fmt=fmt,type=type,data=data,offs=offs,npiece=np,cf=rf)
  +
  +      if(E_IO == 0) then
  +        allocate(var(nx1:nx2,ny1:ny2,nz1:nz2), 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)
  +            if (trim(adjustlt(Upper_Case(fmt)))/='BINARY' .or. &
  +                trim(adjustlt(Upper_Case(type)))/='INT16') then
  +              E_IO = -1_I4P
  +            else
  +              data=trim(adjustlt(data))
  +              allocate(dI1P((nx2-nx1+1)*(ny2-ny1+1)*(nz2-nz1+1)*int(BYI2P,I4P)+int(BYI4P,I4P)))
  +              call b64_decode(code=data,n=dI1P); if(allocated(data)) deallocate(data)
  +              ! Unpack data [1xI4P,(nx2-nx1+1)*(ny2-ny1+1)*(z2-nz1+1)xI2P]
  +              N_byte =  transfer(dI1P(1:int(BYI4P,I4P)),N_byte)
  +              s = size(transfer(dI1P(int(BYI4P,I4P)+1:),XYZp)); allocate(XYZp(1:s))
  +              XYZp = transfer(dI1P(int(BYI4P,I4P)+1:),XYZp); if(allocated(dI1P)) deallocate(dI1P)
  +              var(:,:,:) = reshape(XYZp, (/nx2-nx1+1,ny2-ny1+1,nz2-nz1+1/))
  +              if(allocated(XYZp)) deallocate(XYZp)
  +            endif
  +
  +          case(raw)
  +            if(E_IO == 0) then
  +              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
         endif
  -
  -    endselect
  -  endif
  +  endselect
     !---------------------------------------------------------------------------------------------------------------------------------
     end function
   
  
positive.code:   character(len=:), allocatable          :: fmt
    character(len=:), allocatable          :: type
    character(len=:), allocatable          :: data
    integer(I4P)                           :: np, offs
    integer(I4P)                           :: N_Byte
    !---------------------------------------------------------------------------------------------------------------------------------
  
    !---------------------------------------------------------------------------------------------------------------------------------
  
negative.code:       E_IO = move(inside='UnstructuredGrid', to_find='Piece', repeat=np, buffer=s_buffer)
        inquire(unit=vtk(rf)%u, pos=pos, iostat=E_IO) !annotate the current position in the file
        call get_int(buffer=s_buffer, attrib='NumberOfCells', val=NC, E_IO=E_IO)
        ! get ascii array offsets
        allocate(offset(NC), stat=E_IO)
        E_IO = search(from=pos,inside='Cells', to_find='DataArray', with_attribute='Name', of_value='Offsets', &
                      buffer=s_buffer, content=data)
  