anchor.diff:         end if
       end do
   
  -  end subroutine compute_vxz
  -
  -  subroutine div_vxz(n_frames)
  -    integer, intent(in) :: n_frames
  +  end subroutine compute_vxz_and_vxyz
  +
  +  subroutine div_vxz()
   
       integer :: i, j
       do i = 1, L(1)
          do j = 1, L(3)
             if (v_xz_count(j, i) > 0) then
  -             v_xz(:, j, i) = v_xz(:, j, i) / (v_xz_count(j, i)*n_frames)
  +             v_xz(:, j, i) = v_xz(:, j, i) / v_xz_count(j, i)
             end if
          end do
       end do
   
     end subroutine div_vxz
   
  +  subroutine div_vxyz()
  +
  +    integer :: i, j
  +    do i = 1, L(1)
  +       do j = 1, L(2)
  +          do k = 1, L(3)
  +             if (v_xyz_count(k, j, i) > 0) then
  +                v_xyz(:, k, j, i) = v_xyz(:, k, j, i) / v_xyz_count(k, j, i)
  +             end if
  +          end do
  +       end do
  +    end do
  +
  +  end subroutine div_vxyz
  +
   end program setup_sphere_thermo_trap
  
positive.code: 
    end subroutine rescale_at_walls
  
    subroutine compute_vxz
      integer :: i, s, ix, iy, iz
  
      do i = 1, solvent%Nmax
  
negative.code:        ix = modulo(floor(solvent%pos(1,i)/solvent_cells%a), L(1)) + 1
         iy = modulo(floor(solvent%pos(2,i)/solvent_cells%a), L(2)) + 1
         iz = modulo(floor(solvent%pos(3,i)/solvent_cells%a), L(3)) + 1
         v_xyz_count(iz, iy, ix) = v_xyz_count(iz, iy, ix) + 1
         v_xyz(:, iz, iy, ix) = v_xyz(:, iz, iy, ix) + solvent%vel(:, i)
         if ( (iy == L(2)/2) .or. (iy == 1+L(2)/2) ) then
            v_xz_count(iz, ix) = v_xz_count(iz, ix) + 1
            v_xz(1, iz, ix) = v_xz(1, iz, ix) + solvent%vel(1, i)
  