anchor.diff:              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 )
  +              if( .not. associated(surf(i)%neighbor) ) then
  +                allocate( surf(i)%neighbor(N_NEIGHBOR_MAX_INIT) )
  +                surf(i)%n_neighbor_max = N_NEIGHBOR_MAX_INIT
  +                surf(i)%neighbor(1) = j
  +              else if( surf(i)%n_neighbor > surf(i)%n_neighbor_max ) then
  +                oldsize = surf(i)%n_neighbor_max
  +                newsize = oldsize * 2
  +                dumarray => surf(i)%neighbor
  +                allocate( surf(i)%neighbor(newsize) )
  +                surf(i)%n_neighbor_max = newsize
  +                do k=1,oldsize
  +                  surf(i)%neighbor(k) = dumarray(k)
  +                enddo
  +                surf(i)%neighbor(oldsize+1) = j
  +                deallocate( dumarray )
  +              else
  +                surf(i)%neighbor(surf(i)%n_neighbor) = j
  +              endif
   
               endif
             enddo
           enddo
         enddo
       enddo
  +!$omp end parallel do
  +
     end subroutine
   
     !> Tracing next contact position
  
positive.code: integer, parameter       :: l_max_surface_node =20
  integer, parameter       :: l_max_elem_node = 100
  
      !> 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, pointer     :: neighbor(:)=>null()  !< index(global) of neighborhood tSurfElement
      end type
  
  
negative.code:     surf%nodes(1:n)=nodes(1:n)
      n=getNumberOfSubface( outtype )
      surf%n_neighbor = 0
      surf%n_neighbor_max = 0
    end subroutine
  
    !> Memeory management subroutine
  