anchor.diff:      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))
  +!  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))
  +!  print*,"before collision", maxval(cvx), cpx(maxloc(cvx)), cpy(maxloc(cvx)),cpz(maxloc(cvx))
   
     ! Resolve collisions
     do i = 1, n
  
positive.code:   do subiter = 1, ITERATIONS_PER_STEP
      ! Compute lamnda
      call compute_lambda
    print*,'after lambda', maxval(cvx), maxval(cvy), maxval(cvz)
      ! Jiggle the particles
      call particle_jiggle
    print*,'after jiggle', maxval(cvx), maxval(cvy), maxval(cvz)
    end do
  
    ! Confine vorticity, update vel, and calculate force
    if(USE_VORTICITY_CONFINEMENT) then
  
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
  
  