Computational methodology for exact quantum 6D calculations of the vibrational eigenstates, energy levels, and wave functions of a diatomic molecule adsorbed on a rigid corrugated surface is presented. It is intended for adsorbates executing coupled, strongly anharmonic large-amplitude vibrations. Surface nonrigidity is introduced in an approximate way, by means of a simplified surface-mass model. Using this methodology, we calculate the vibrational levels of CO/Cu(100) for all four isotopomers of CO, (CO)-C-12-O-16, (CO)-C-13-O-16, (CO)-C-12-O-18, and (CO)-C-13-O-18. The empirical potential by Tully and co-workers [J. C. Tully, M. Gomez, and M. Head-Gordon, J. Vac. Sci. Technol. A 11, 1914 (1993)] is employed. Our calculated fundamental frequencies of CO/Cu(100) vibrations are compared to those from earlier theoretical treatments on the same potential, as well as with the experimental frequencies and isotope frequency shifts. In addition to 6D calculations, we perform 5D (rigid CO) and 4D (fixed-site) quantum calculations, which provide information about the couplings among the vibrational modes of CO on Cu(100). Excited levels of the lowest-frequency in-plane (doubly degenerate) frustrated translation mode are analyzed and assigned.