The energy levels and wavefunctions including the two lowest-lying levels, namely ⁴A₂ and _a²E, for which reliable experimental data exist for Cr³⁺ ions at C₃ symmetry sites in LiNbO₃, are calculated using the complete matrix diagonalization method within the 3d³ configuration. The Hamiltonian considered includes the electrostatic term, the Trees correction, the spin-orbit interaction and the crystal-field interaction. The role of the additional low-symmetry crystal-field term B₄^-3O₄^-3 (in the Stevens operator notation), neglected in the C_3v approximation used so far in the literature, is studied. The superposition model is developed for 3d³ ions at C₃ symmetry sites and applied to study the site occupancy of Cr³⁺ in LiNbO₃. Analysis of the optical data indicates that Cr³⁺ ions substitute at Nb sites and Li sites simultaneously. The present considerations offer an improvement over the earlier approximations using C_3v symmetry only. The zero-field splitting predicted by the crystal-field calculations for Cr³⁺ at the Nb site matches the experimental value from EPR studies very well. This is contrary to the earlier prediction by the superposition model analysis of the spin-Hamiltonian parameters indicating that the zero-field splitting for Cr³⁺ ions at Li sites matches the experimental zero-field splitting better than that for Cr³⁺ at Nb sites. Since the present calculations involve fitting not only the zero-field splitting but also the energies of the _a²E state, the present predictions may be more reliable than the previous predictions.