anchor.diff:    return
   end subroutine compute_visc
   
  +subroutine compute_SF
  +  use posvel
  +  use forces
  +!http://www.idi.ntnu.no/~elster/master-studs/fossum/fossum-fall2010-proj.pdf
  +
  +
  +
  +  return
  +end subroutine compute_SF
  +
   subroutine ext_forces
     use forces
     use constants
  
positive.code:   integer :: binx, biny, binz
    integer :: cbinx, cbiny, cbinz
    integer :: stx, sty, stz
    real(WP) :: distx, disty, distz, tdist
  
    part_count = 0
    binpart = 0
  
negative.code:         do stz = -1, 1
            cbinz = binz + stz
            if(cbinz.lt.1 .or. cbinz.gt.nbinz) cycle
            do j = 1, part_count(cbinx, cbiny, cbinz)
              p = binpart(j,cbinx,cbiny,cbinz)
              distx = px(i) - px(p)
              disty = py(i) - py(p)
              distz = pz(i) - pz(p)
              tdist2 = distx*distx + disty*disty + distz*distz
              if(tdist2 .lt. h2 .and. p.ne.i) then
                q = q + 1
                nbs(i,q) = binpart(j,cbinx,cbiny,cbinz)
              end if
  