anchor.diff:      l = nc(i)
       do j = 1, l
         nb = nbs(i,j)
  -
  +      if(nb.lt.i) cycle
         dx = px(i) - px(nb)
         dy = py(i) - py(nb)
         dz = pz(i) - pz(nb)
  +      dvx = vx(i) - vx(nb)
  +      dvy = vy(i) - vy(nb)
  +      dvz = vz(i) - vz(nb)
  +
         r2 = dx*dx + dy*dy + dz*dz
         m = sqrt(r2/h2)
         q = 1.0_WP - m
   
   
  -      k = c / rho(nb) * q &
  +      kpres = cpres / rho(nb) * q &
             * (rho(i)+rho(nb)-2.0_WP*PARTICLE_DENSITY) &
             * q / m
   
  -      sx = sx + k*dx
  -      sy = sy + k*dy
  -      sz = sz + k*dz
  +      kvisc = cvisc / rho(nb) * q
  +
  +
  +      sx = sx + kpres*dx+kvisc*dvx
  +      fx(nb) = fx(nb) - kpres*dx-kvisc*dvx
  +      sy = sy + kpres*dy+kvisc*dvy
  +      fy(nb) = fy(nb) - kpres*dy-kvisc*dvy
  +      sz = sz + kpres*dz+kvisc*dvz
  +      fz(nb) = fz(nb) - kpres*dz-kvisc*dvz
   
       end do
  -
  -    fx(i) = sx/rho(i)
  -    fy(i) = sy/rho(i)
  -    fz(i) = sz/rho(i)
  -
  -  end do
  -
  -  return
  -end subroutine compute_pressure
  -
  -subroutine compute_visc
  -  use constants
  -  use neighbors
  -  use posvel
  -  use forces
  -  use state
  -  implicit none
  -  integer :: i, j, l, nb
  -  real(WP) :: sx, sy, sz
  -  real(WP) :: dvx, dvy, dvz
  -  real(WP) :: dx, dy, dz
  -  real(WP) :: k, c, r2
  -
  -  c = -45.0_WP*mass/PI/h5*VISCOSITY
  -
  -  do i = 1, n
  -    sx = 0.0_WP
  -    sy = 0.0_WP
  -    sz = 0.0_WP
  -
  -    ! Compute viscous
  -    l = nc(i)
  -    do j = 1, l
  -      nb = nbs(i,j)
  -
  -      dvx = vx(i) - vx(nb)
  -      dvy = vy(i) - vy(nb)
  -      dvz = vz(i) - vz(nb)
  -      dx = px(i) - px(nb)
  -      dy = py(i) - py(nb)
  -      dz = pz(i) - pz(nb)
  -      r2 = dx*dx + dy*dy + dz*dz
  -
  -      k = c / rho(nb) * (1.0_WP-sqrt(r2/h2))
  -
  -      sx = sx + k*dvx
  -      sy = sy + k*dvy
  -      sz = sz + k*dvz
  -
  -    end do
  -
  -    fx(i) = fx(i) + sx/rho(i)
  -    fy(i) = fy(i) + sy/rho(i)
  -    fz(i) = fz(i) + sz/rho(i)
  -
  -  end do
  -
  -  return
  -end subroutine compute_visc
  +    fx(i) = fx(i)+sx
  +    fy(i) = fy(i)+sy
  +    fz(i) = fz(i)+sz
  +
  +
  +  end do
  +
  +  return
  +end subroutine compute_forces
   
   subroutine compute_SF
     use posvel
  
positive.code:   return
  end subroutine compute_density
  
  subroutine compute_pressure
    use constants
    use neighbors
    use posvel
  
negative.code:   use precision
    implicit none
  
    real(WP) :: neigh_t, dens_t, force_t, ef_t, ps_t, tot_t
  
  end module timers
  
  