anchor.diff:      sub_name = 'define_data_geometry'
       call enter_exit(sub_name,1)
       
  +    if(index(datafile, ".ipdata")> 0) then !full filename is given
  +       readfile = datafile
  +    else ! need to append the correct filename extension
  +       readfile = trim(datafile)//'.ipdata'
  +    endif    !readfile = trim(datafile)//'.ipdata'
  +    open(10, file=readfile, status='old')
  +    
       !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
  
positive.code:     sub_name = 'make_data_grid'
      call enter_exit(sub_name,1)
      
      call triangles_from_surface(num_triangles,num_vertices,surface_elems, &
           triangle,vertex_xyz)
  
      if(offset.gt.0.0_dp)then
  !!! generate within a scaled mesh, then return to original size afterwards
  !!! this option is used when we don't want the seed points too close to the boundary
         scale_mesh = 1.0_dp-(offset/100.0_dp)
         cofm_surfaces = sum(vertex_xyz,dim=2)/size(vertex_xyz,dim=2)
         translate = cofm_surfaces * (scale_mesh - 1.0_dp)
         forall (i=1:num_vertices) vertex_xyz(1:3,i) = &
              vertex_xyz(1:3,i)*scale_mesh + translate(1:3)
      endif
      
  !!! find the bounding coordinates for the surface mesh
      
      min_bound = minval(vertex_xyz,2)
  
negative.code:     if(allocated(data_weight)) deallocate(data_weight)
      allocate(data_weight(3,num_data))
      data_weight(:,1:num_data) = 1.0_dp
  
      deallocate(triangle)
      deallocate(vertex_xyz)
      deallocate(elem_list)
     
      call enter_exit(sub_name,2)
      
    end subroutine make_data_grid
  