anchor.diff:    character(len=:), allocatable       :: fmt
     character(len=:), allocatable       :: type
     character(len=:), allocatable       :: data
  -  integer(I4P)                        :: np, i, j, k, offs, N_Byte, pos
  -
  +  integer(I4P)                        :: np, i, j, k, offs, N_Byte, pos, s
  +  integer(I1P), allocatable           :: dI1P(:)
  +  real(R4P), allocatable              :: XYZp(:)
     !---------------------------------------------------------------------------------------------------------------------------------
   
     !---------------------------------------------------------------------------------------------------------------------------------
  
positive.code:       endif
      end select
    !---------------------------------------------------------------------------------------------------------------------------------
    end function VTK_GEO_XML_STRG_1DAP_R4_READ
  
  
    function VTK_GEO_XML_STRG_1DAP_R8_READ(nx1,nx2,ny1,ny2,nz1,nz2,NN,XYZ,npiece,cf) result(E_IO)
    !---------------------------------------------------------------------------------------------------------------------------------
    !< Function for reading mesh with \b UnstructuredGrid topology (R8P).
    !---------------------------------------------------------------------------------------------------------------------------------
    integer(I4P),           intent(OUT) :: nx1      !< Initial node of x axis.
    integer(I4P),           intent(OUT) :: nx2      !< Final node of x axis.
    integer(I4P),           intent(OUT) :: ny1      !< Initial node of y axis.
    integer(I4P),           intent(OUT) :: ny2      !< Final node of y axis.
    integer(I4P),           intent(OUT) :: nz1      !< Initial node of z axis.
    integer(I4P),           intent(OUT) :: nz2      !< Final node of z axis.
    integer(I4P),           intent(OUT) :: NN       !< Number of nodes
    real(R8P), allocatable, intent(OUT) :: XYZ(:,:)     !< x coordinates
    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        :: aux      !< Auxiliary string.
    character(len=:), allocatable       :: fmt
    character(len=:), allocatable       :: type
    character(len=:), allocatable       :: data
    integer(I4P)                        :: np, i, j, offs, N_Byte, pos
  
    !---------------------------------------------------------------------------------------------------------------------------------
  
    !---------------------------------------------------------------------------------------------------------------------------------
    rf = f
    if (present(cf)) then
      rf = cf ; f = cf
    endif
    np = 1_I4P; if (present(npiece)) np = npiece
    E_IO = -1_I4P
    select case(vtk(rf)%f)         
  
      case(ascii)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside='StructuredGrid', to_find='Piece', repeat=np, cf=rf, buffer=s_buffer) ! find the 'np' piece
        inquire(unit=vtk(rf)%u, pos=pos, iostat=E_IO) !annotate the current position in the file
        call get_char(buffer=s_buffer, attrib='Extent', val=aux, E_IO=E_IO)
        if(E_IO == 0) then      
          read(aux,*, iostat=E_IO) nx1,nx2,ny1,ny2,nz1,nz2
          if(E_IO == 0) then      
            NN = (nx2-nx1+1)*(ny2-ny1+1)*(nz2-nz1+1)
            E_IO = search(from=pos, inside='Points', to_find='DataArray', with_attribute='Name', of_value='Points', &
                          buffer=s_buffer,content=data)
            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
              allocate(XYZ(3,NN), stat=E_IO)  
              read(data, fmt=*, iostat=E_IO) ((XYZ(i,j),i=1,3),j=1,NN) !get ascii array
            endif
            if(allocated(data)) deallocate(data)
          endif
        endif
  
      case(binary)
  ! not implemented
  
      case(raw)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside='StructuredGrid', to_find='Piece', repeat=np, cf=rf, buffer=s_buffer) ! find the 'np' piece
        inquire(unit=vtk(rf)%u, pos=pos, iostat=E_IO) !annotate the current position in the file
        call get_char(buffer=s_buffer, attrib='Extent', val=aux, E_IO=E_IO)
        if(E_IO == 0) then      
          read(aux,*, iostat=E_IO) nx1,nx2,ny1,ny2,nz1,nz2
          if(E_IO == 0) then      
            NN = (nx2-nx1+1)*(ny2-ny1+1)*(nz2-nz1+1)
            E_IO = search(from=pos, inside='Points', to_find='DataArray', with_attribute='Name', of_value='Points', &
                          buffer=s_buffer)
            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)))/='APPENDED' .or. &
                trim(adjustlt(Upper_Case(type)))/='FLOAT64') then
              E_IO = -1_I4P 
            else
              allocate(XYZ(3,NN), stat=E_IO)        
              read(unit=vtk(rf)%u, iostat=E_IO, pos = vtk(rf)%ioffset+offs) N_Byte, ((XYZ(i,j),i=1,3),j=1,NN) !get appended array
            endif
          endif
        endif
      end select
    !---------------------------------------------------------------------------------------------------------------------------------
    end function VTK_GEO_XML_STRG_1DAP_R8_READ
  
  
  
negative.code:       endif
  
      case(binary)
  
        rewind(unit=vtk(rf)%u, iostat=E_IO)
        E_IO = move(inside='StructuredGrid', to_find='Piece', repeat=np, cf=rf, buffer=s_buffer) ! find the 'np' piece
        inquire(unit=vtk(rf)%u, pos=pos, iostat=E_IO) !annotate the current position in the file
        call get_char(buffer=s_buffer, attrib='Extent', val=aux, E_IO=E_IO)
        if(E_IO == 0) then      
          read(aux,*, iostat=E_IO) nx1,nx2,ny1,ny2,nz1,nz2
          if(E_IO == 0) then      
            NN = (nx2-nx1+1)*(ny2-ny1+1)*(nz2-nz1+1)
            E_IO = search(from=pos, inside='Points', to_find='DataArray', with_attribute='Name', of_value='Points', &
                          buffer=s_buffer,content=data)
            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)))/='BINARY' .or. &
                trim(adjustlt(Upper_Case(type)))/='FLOAT32') then
              E_IO = -1_I4P 
            else
              ! Decode base64 packed data
              data=trim(adjustlt(data))
              allocate(dI1P(3*NN*int(BYR4P,I4P)+int(BYI4P,I4P)))
              call b64_decode(code=data,n=dI1P); if(allocated(data)) deallocate(data)
              ! Unpack data [1xI4P,3*NNxR4P]
              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)
              allocate(X(NN), Y(NN), Z(NN), stat=E_IO)
              do i=1,NN; X(i)=XYZp(i*3-2); Y(i)=XYZp(i*3-1); Z(i)=XYZp(i*3); enddo
              if(allocated(XYZp)) deallocate(XYZp)
            endif
          endif
        endif
  
      case(raw)
  
  