A new closure to the Ornstein-Zernike (OZ) equation is proposed for ionic liquids and is investigated for primitive models of high valency (2:2) aqueous electrolyte solutions and molten salts. The new closure, which is related to an earlier closure for the soft-sphere case proposed by Ichiye and Haymet, may be viewed as a prescription for the so-called "bridge functions." These functions are approximated by zero in the hypernetted-chain (HNC) closure which is generally used for ionic systems. In both the new closure and the soft-sphere closure, the recognition that the unlike bridge function is opposite in sign from the like bridge function leads to an approximation for these missing graphs by adding (for the unlike case) or subtracting (for the like case) a set of graphs similar to those used in Percus-Yevick theory to the HNC equation. Compared to the HNC closure, the pair correlation functions predicted for primitive models by the new closure are generally in much better agreement with Monte Carlo (MC) simulations of molten salts and aqueous 2:2 electrolytes. The fundamental improvement of this paper over the Ichiye-Haymet work is that the separation of long- and short-range part of c(r) for the hard-sphere case is clearly defined, whereas it was done numerically for the soft-sphere case. Moreover, the present theory is in better agreement with MC simulations both in the molten salt as well as in the dilute solution regimes than the soft-sphere case. Finally, a study was made of the transition of the like charge pair correlation functions from monotonic behavior at low densities to a nonmonotonic behavior at high densities. The new closure clearly predicts such a transition region at concentrations near 0.02 M and temperatures near 314 K. There is also a region below 0.02 M and 314 K where the new closure fails to converge. Compared to MC simulations, the critical region predicted by the new closure appears to be a lower estimate. However, for the HNC closure there is only a remote possibility of such a transition region since the correlation functions are nonmonotonic even at lower concentrations, a feature which is corrected in the new theory.