anchor.diff:             ijunk  =intfromfile(6)
              lcsyst =intfromfile(7)
              allocate (ientp(neltp,nshl))
  -c           read(igeomBAK) ientp
              iientpsiz=neltp*nshl
   
              if (neltp==0) then
  
positive.code:         use pointer_data
          use phio
          use iso_c_binding
  c
          include "common.h"
  !MR CHANGE
          include "mpif.h" !Required to determine the max for itpblk
  !MR CHANGE END
  c
          integer, target, allocatable :: ientp(:,:)
          integer mater(ibksz)
          integer, target :: intfromfile(50) ! integers read from headers
          character*255 fname1
  
  cccccccccccccc New Phasta IO starts here ccccccccccccccccccccccccc
  
          integer :: descriptor, descriptorG, GPID, color, nfiles
          integer ::  numparts, writeLock
          integer :: ierr_io, numprocs
  !MR CHANGE
          integer, target :: itpblktot,ierr,iseven
  !MR CHANGE END
          character*255 fname2
   
          character(len=30) :: dataInt
          dataInt = c_char_'integer'//c_null_char
  
  !THIS NEEDS TO BE CLEANED - MR
          nfiles = nsynciofiles
  !        nfields = nsynciofieldsreadgeombc
          numparts = numpe !This is the common settings. Beware if you try to compute several parts per process
  
  
          ione=1
          itwo=2
          iseven=7
          ieleven=11
  
  cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  
  c
          iel=1
          itpblk=nelblk
  !MR CHANGE
  
          ! Get the total number of different interior topologies in the whole domain. 
          ! Try to read from a field. If the field does not exist, scan the geombc file.
  
negative.code:         call phio_readheader(fhandle,
       &   c_char_'total number of interior tpblocks' // char(0),
       &   c_loc(itpblktot), ione, dataInt, iotype) 
  
          if (itpblktot == -1) then 
            ! The field 'total number of different interior tpblocks' was not found in the geombc file.
  