anchor.diff:      call compute_omega
       call vorticity_force
     end if
  -!  print*,'before visc', maxval(cvx), maxval(cvy), maxval(cvz)
  +  print*,'before visc', maxval(cvx), maxval(cvy), maxval(cvz)
     ! Update particle positions due to vorticity
     call viscosity_update
  -!  print*,'update with visc', maxval(cvx), maxval(cvy), maxval(cvz)
  +  print*,'update with visc', maxval(cvx), maxval(cvy), maxval(cvz)
   
     ! Update particle position
     call pos_update
  
positive.code:     dpy(i) = sy / PARTICLE_DENSITY
      dpz(i) = sz / PARTICLE_DENSITY
    end do
    print*,"before update", maxval(cvx), cpx(maxloc(cvx)), cpy(maxloc(cvx)),cpz(maxloc(cvx))
      ! Update candidate position
      cpx(:) = cpx(:) + dpx(:)
      cpy(:) = cpy(:) + dpy(:)
      cpz(:) = cpz(:) + dpz(:)
  
      ! Update candidate velocity
      cvx(:) = (cpx(:) - px(:)) / dt
      cvy(:) = (cpy(:) - py(:)) / dt
      cvz(:) = (cpz(:) - pz(:)) / dt
  
    print*,"before collision", maxval(cvx), cpx(maxloc(cvx)), cpy(maxloc(cvx)),cpz(maxloc(cvx))
  
    ! Resolve collisions
    do i = 1, n
  
negative.code: 
    real(WP), parameter :: ARTIFICIAL_VISCOSITY = 1.0e-2_WP
    real(WP), parameter :: VORTICITY_COEFFICIENT = 1.0e-4_WP
    logical,  parameter :: USE_VORTICITY_CONFINEMENT = .true.
  
    real(WP) :: PRESSURE_RADIUS_FACTOR
  
  