anchor.diff:      surf%nodes(1:n)=nodes(1:n)
       n=getNumberOfSubface( outtype )
       surf%n_neighbor = 0
  +    surf%n_neighbor_max = 0
     end subroutine
   
     !> Memeory management subroutine
  
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:   !> 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)
  