anchor.diff:          call getShapeFunc( etype, opos(:), shapefunc )
           do j=1,nn
             iSS = contact%master(omaster)%nodes(j)
  -          B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)-nrlforce*shapefunc(j)*odirec
  +          ! B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)-nrlforce*shapefunc(j)*odirec
  +          idx0 = 3*(iSS-1)
  +          do i=1,3
  +!$omp atomic
  +            B(idx0+i) = B(idx0+i)-nrlforce*shapefunc(j)*odirec(i)
  +          enddo
           enddo
           if( fcoeff/=0.d0 ) then
             do j=1,nn
  
positive.code:       real(kind=kreal)    :: dum, dxi(2), shapefunc(l_max_surface_node)
        real(kind=kreal)    :: metric(2,2), dispmat(2,l_max_elem_node*3+3)
        real(kind=kreal)    :: fric(2), f3(l_max_elem_node*3+3)
  
        slave = contact%slave(lslave)
        fcoeff = contact%fcoeff
  
negative.code:                        " along direction ", contact%states(i)%direction
              exit
            enddo
            deallocate(indexMaster)
          endif
        enddo
  !$omp end parallel do
  
    end subroutine scan_contact_state
  
    !> Calculate averaged nodal normal
  