  (ABRIDGED) We present a new Schwarzschild orbit-superposition code designed
to model discrete datasets composed of velocities of individual kinematic
tracers in a dynamical system. This constitutes an extension of previous
implementations that can only address continuous data in the form of (the
moments of) velocity distributions, thus avoiding potentially important losses
of information due to data binning. Furthermore, the code can handle any
combination of available velocity components, i.e., only line-of-sight
velocities, only proper motions, or a combination of both. It can also handle a
combination of discrete and continuous data. The code finds the distribution
function (DF, a function of the three integrals of motion E, Lz, and I3) that
best reproduces the available kinematic and photometric observations in a given
axisymmetric gravitational potential. The fully numerical approach ensures
considerable freedom on the form of the DF f(E,Lz,I3). This allows a very
general modeling of the orbital structure, thus avoiding restrictive
assumptions about the degree of (an)isotropy of the orbits. We describe the
implementation of the discrete code and present a series of tests of its
performance based on the modeling of simulated datasets generated from a known
DF. We find that the discrete Schwarzschild code recovers the original orbital
structure, M/L ratios, and inclination of the input datasets to satisfactory
accuracy, as quantified by various statistics. The code will be valuable, e.g.,
for modeling stellar motions in Galactic globular clusters, and those of
individual stars, planetary nebulae, or globular clusters in nearby galaxies.
This can shed new light on the total mass distributions of these systems, with
central black holes and dark matter halos being of particular interest.
