An algorithm is developed for passivity preserving model reduction of LTI systems. The derivation is justified analytically and implementation schemes are developed for both medium scale (dense) and large scale (sparse) applications. The algorithm is based upon interpolation of specified spectral zeros of the original transfer function to produce a reduced transfer function that has the specified roots as its spectral zeros. These interpolation conditions are satisfied through the computation of a basis for a selected invariant subspace of a certain blocked matrix which has the spectral zeros as its spectrum. It is shown that this procedure indirectly solves the associated controllability and observability Riccati equations and how to select the interpolation points to give maximal or minimal solutions of these equations. From these, a balancing transformation may be constructed to give a reduced model that is balanced as well as passive and stable.