anchor.diff:  
   !!!#############################################################################
     
  -  subroutine make_data_grid(surface_elems,spacing,to_export,filename,groupname)
  +  subroutine make_data_grid(surface_elems, offset, spacing, filename, groupname)
       !*make_data_grid:* makes a regularly-spaced 3D grid of data points to
       ! fill a bounding surface 
       !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_MAKE_DATA_GRID" :: MAKE_DATA_GRID
  +
  +    use exports,only: export_triangle_elements,export_triangle_nodes
       
       integer,intent(in) :: surface_elems(:)
  -    real(dp),intent(in) :: spacing
  -    logical,intent(in) :: to_export
  +     real(dp),intent(in) :: offset, spacing
  +    logical :: to_export = .true.
       character(len=*),intent(in) :: filename
       character(len=*),intent(in) :: groupname
       ! Local Variables
       integer :: i,j,k,ne,nj,nline,nn,num_data_estimate,num_triangles,num_vertices
  -    integer,allocatable :: triangle(:,:)
  -    real(dp) :: cofm_surfaces(3),boxrange(3),max_bound(3),min_bound(3), &
  -         offset=-2.0_dp,point_xyz(3),scale_mesh,translate(3)
  +    integer,allocatable :: elem_list(:),triangle(:,:)
  +    real(dp) :: cofm1(3),cofm2(3),boxrange(3),max_bound(3),min_bound(3), &
  +         point_xyz(3),scale_mesh
       real(dp),allocatable :: data_temp(:,:),vertex_xyz(:,:)
       logical :: internal
       character(len=1) :: char1
  
positive.code:     if(allocated(data_weight)) deallocate(data_weight)
      allocate(data_weight(3,num_data))
      data_weight(:,1:num_data) = 1.0_dp
      if(offset .gt. 0.0_dp)then
  !!! return the triangles to their original size and location
         forall (i=1:3) vertex_xyz(1:3,i) = (vertex_xyz(1:3,i) - &
              translate(1:3))/scale_mesh
      endif
      
      if(to_export)then
  !!! export vertices as nodes
         writefile = trim(filename)//'.exnode'
         open(10, file = writefile, status='replace')
         !**    write the group name
         write(10,'( '' Group name: '',A)') trim(groupname)
         !*** Exporting Geometry
         !*** Write the field information
         write(10,'( '' #Fields=1'' )')
         write(10,'('' 1) coordinates, coordinate, rectangular cartesian, '',&
              &''#Components=3'')')
         do nj=1,3
            if(nj.eq.1) write(10,'(2X,''x.  '')',advance="no")
            if(nj.eq.2) write(10,'(2X,''y.  '')',advance="no")
            if(nj.eq.3) write(10,'(2X,''z.  '')',advance="no")
            write(10,'(''Value index='',I2,'', #Derivatives='',I1)', &
                 advance="no") nj,0
            write(10,'()')
         enddo
         do i = 1,num_vertices
            !***    write the node
            write(10,'(1X,''Node: '',I12)') i
            write(10,'(2x,3(f12.6))') vertex_xyz(:,i)
         enddo
         close(10)
         
  !!! export the triangles as surface elements
         writefile = trim(filename)//'.exelem'
         open(10, file=writefile, status='replace')
         !**     write the group name
         write(10,'( '' Group name: '',a10)') groupname
         !**     write the lines
         write(10,'( '' Shape. Dimension=1'' )')
         nline=0
         do ne = 1,num_triangles
            write(10,'( '' Element: 0 0 '',I5)') nline+1
            write(10,'( '' Element: 0 0 '',I5)') nline+2
            write(10,'( '' Element: 0 0 '',I5)') nline+3
            nline=nline+3
         enddo !ne
         
         !**        write the elements
         write(10,'( '' Shape. Dimension=2, line*line'' )')
         write(10,'( '' #Scale factor sets=1'' )')
         write(10,'( '' l.Lagrange*l.Lagrange, #Scale factors=4'' )')
         write(10,'( '' #Nodes= '',I2 )') 4
         write(10,'( '' #Fields=1'' )')
         write(10,'( '' 1) coordinates, coordinate, rectangular cartesian, #Components=3'')')
         
         do nj=1,3
            if(nj==1) char1='x'; if(nj==2) char1='y'; if(nj==3) char1='z';
            write(10,'(''  '',A2,''. l.Lagrange*l.Lagrange, no modify, standard node based.'')') char1
            write(10,'( ''   #Nodes=4'')')
            do nn=1,4
               write(10,'(''   '',I1,''. #Values=1'')') nn
               write(10,'(''     Value indices: '',I4)') 1
               write(10,'(''     Scale factor indices: '',I4)') nn
            enddo !nn
         enddo !nj
         
         nline=0
         do ne=1,num_triangles
            !**         write the element
            write(10,'(1X,''Element: '',I12,'' 0 0'' )') ne
            !**          write the faces
            WRITE(10,'(3X,''Faces: '' )')
            WRITE(10,'(5X,''0 0'',I6)') nline+1
            WRITE(10,'(5X,''0 0'',I6)') nline+2
            WRITE(10,'(5X,''0 0'',I6)') nline+3
            WRITE(10,'(5X,''0 0'',I6)') 0
            nline=nline+3
            !**          write the nodes
            write(10,'(3X,''Nodes:'' )')
            write(10,'(4X,4(1X,I12))') triangle(:,ne),triangle(3,ne)
            !**          write the scale factors
            write(10,'(3X,''Scale factors:'' )')
            write(10,'(4X,4(1X,E12.5))') 1.0_dp,1.0_dp,1.0_dp,1.0_dp
         enddo
         close(10)
      endif
      
      deallocate(triangle)
      deallocate(vertex_xyz)
      
      call enter_exit(sub_name,2)
      
    end subroutine make_data_grid
  
negative.code:     end do read_number_of_elements
      
      num_elems_2d=number_of_elements
      if(allocated(elems_2d))then
         deallocate(elems_2d)
         deallocate(elem_nodes_2d)
         deallocate(elem_versn_2d)
      endif
      allocate(elems_2d(num_elems_2d))
      allocate(elem_nodes_2d(4,num_elems_2d))
      allocate(elem_versn_2d(4,num_elems_2d))
      
      ne = 0
      
  