anchor.diff:      elem_cnct(-1:1,0:2,0:num_elems)=enodes_temp2(-1:1,0:2,0:num_elems)
       deallocate(enodes_temp2)
       elem_cnct(-1:1,0:2,num_elems+1:num_elems_new) = 0
  -
  +    
       allocate(enodes_temp(num_ord,num_elems))
       enodes_temp=elem_ordrs
       deallocate(elem_ordrs)
  
positive.code:     if(.not.allocated(elem_field)) allocate(elem_field(num_ne,num_elems))
      if(.not.allocated(elem_direction)) allocate(elem_direction(3,num_elems))
      if(.not.allocated(expansile)) allocate(expansile(num_elems))
  
  !!! initialise array values
      nodes = 0 !initialise node index values
      node_xyz = 0.0_dp !initialise
  
negative.code:     end do read_a_node
      
      close(10)
      
      call enter_exit(sub_name,2)
      
    end subroutine define_node_geometry_2d
  
  !!!#############################################################################
  
    subroutine define_data_geometry(datafile)
      !*define_data_geometry:* reads data points from a file
      !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_DEFINE_DATA_GEOMETRY" :: DEFINE_DATA_GEOMETRY
  
      character(len=*) :: datafile
      ! Local variables
      integer :: iend,ierror,length_string,ncount,nj,itemp
      character(len=132) :: buffer,readfile
      character(len=60) :: sub_name
  
      ! --------------------------------------------------------------------------
  
      sub_name = 'define_data_geometry'
      call enter_exit(sub_name,1)
      
      !set the counted number of data points to zero
      ncount = 0
      
      !readfile = trim(datafile)//'.ipdata'
      open(10, file=datafile, status='old')
      read(unit=10, fmt="(a)", iostat=ierror) buffer
      
  !!! first run through to count the number of data points
      read_line_to_count : do
         read(unit=10, fmt="(a)", iostat=ierror) buffer
  