The aim of this work is twofold. In the first part, a detailed description is given of a specific vibrational model, designed for calculations on the vibrational energy levels in benzene and benzene isotopic species of D-6h, symmetry. For the description of the C-H stretch system in benzene, a local mode (LM) formalism was applied, while, for the remaining non-C-H stretch vibrations, a symmetrized mode (SM) treatment was applied: this was called the combined LM/SM model. The model is based on a set of complex symmetrized curvilinear vibrational coordinates, which can be expressed as simple linear combinations of Whiffen's coordinates. The description in terms of complex symmetrized coordinates and wave functions allows for the construction of a separable symmetrized infinite-dimensional vibrational basis set, which is of crucial importance for large-scale calculations. In the second part of this work, using the described complex symmetrized LM/SM vibrational model, calculations have been carried out on a large number of vibrational energy levels of four benzene D-6h, isotopomers. The aim of the calculations was to redetermine a reliable set of harmonic force constants for benzene. Some of the force constant values obtained in the present work are substantially different from previous determinations by other authors. Using the presently determined set of harmonic force constant values in the calculations, a very good fit has been obtained to a large number of experimentally measured vibrational (both fundamental and overtone) energy levels of various symmetries, belonging to all four D-6h benzene isotopomers: C6H6, C6D6, (C6H6)-C-13, (C6D6)-C-13.