anchor.diff:  integer, parameter       :: l_max_surface_node =20
   integer, parameter       :: l_max_elem_node = 100
   
  +integer, parameter       :: N_NEIGHBOR_MAX_INIT = 8
  +
       !> Structure to define surface group
       type tSurfElement
         integer              :: eid                  !< elemental index(global)
         integer              :: etype                !< type of surface element
         integer, pointer     :: nodes(:)=>null()     !< nodes index(global)
  -      integer              :: n_neighbor           !< size of neighber(:)
  +      integer              :: n_neighbor           !< number of neighborhood tSurfElement
  +      integer              :: n_neighbor_max       !< maximum size of neighber(:)
         integer, pointer     :: neighbor(:)=>null()  !< index(global) of neighborhood tSurfElement
       end type
   
  
positive.code:             if( nd1==nd2 ) then
                surf(i)%n_neighbor = surf(i)%n_neighbor+1
  
        if( .not. associated(surf(i)%neighbor) ) then
          allocate( surf(i)%neighbor(1) )
          surf(i)%neighbor(1) = j
        else
          oldsize = size( surf(i)%neighbor )
          allocate( dumarray(oldsize) )
          do k=1,oldsize
            dumarray(k) = surf(i)%neighbor(k)
          enddo
          deallocate( surf(i)%neighbor )
          allocate( surf(i)%neighbor(oldsize+1) )
          do k=1,oldsize
            surf(i)%neighbor(k) = dumarray(k)
          enddo
          surf(i)%neighbor(oldsize+1) = j
        endif
        if( associated(dumarray) ) deallocate( dumarray )
  
              endif
            enddo
          enddo
        enddo
      enddo
    end subroutine
  
    !> Tracing next contact position
  
negative.code:   !> Find neighboring surface elements
    subroutine find_surface_neighbor( surf )
      use m_utilities
      use hecmw_util
      type(tSurfElement), intent(inout) :: surf(:)
      integer :: i, j, ii,jj, nd1, nd2, nsurf
      integer :: k, oldsize, newsize
      integer, pointer :: dumarray(:) => NULL()
  
      nsurf = size(surf)
  
  !$omp parallel do default(none),private(i,ii,nd1,j,jj,nd2,oldsize,newsize,dumarray,k), &
  !$omp&shared(nsurf,surf)
      do i=1,nsurf
        do ii=1, size(surf(i)%nodes)
          nd1 = surf(i)%nodes(ii)
  