anchor.diff:            * (rho(i)+rho(nb)-2.0_WP*PARTICLE_DENSITY) &
             * q / m
   
  -      sx = sx + (k+scorr)*dx
  -      sy = sy + (k+scorr)*dy
  -      sz = sz + (k+scorr)*dz
  +      sx = sx + k*dx
  +      sy = sy + k*dy
  +      sz = sz + k*dz
   
       end do
   
  
positive.code: subroutine compute_SF
    use posvel
    use forces
  !http://www.idi.ntnu.no/~elster/master-studs/fossum/fossum-fall2010-proj.pdf
  
  
  
    return
  
negative.code: 
  
    c = 45.0_WP*mass/PI/h5*BULK_MODULUS*0.5_WP
  
    do i = 1, n
      sx = 0.0_WP
      sy = 0.0_WP
      sz = 0.0_WP
  
      ! Compute pressure
      l = nc(i)
  