A grid-based method designed to refine adiabatic potential energy surfaces (PES) of molecules via minimizing a suitable objective function is described. The objective function contains deviations from the reference (experimental) (ro)-vibrational energy levels and is based on PES correction values determined at the grid points within a discrete-variable-representation nuclear-motion algorithm and first-order perturbation theory (PT). The proposed PES refinement technique is tested on the ground electronic state of the MgH molecule. The large number of numerical test results obtained suggest the following: (1) first-order PT is able to yield accurate correction values at the grid points representing the PES, and for practical cases there seems to be no need to go to higher orders of PT; (2) with the number of grid points greatly exceeding the number of experimental energy levels included in the refinement procedure, terms additional to the "obs-calc" term, including numerical first and second derivatives of the correction surface, are necessary in the objective function to arrive at a physically meaningful, "smooth" correction surface; (3) for a given J rotational quantum number, the corrected PES is able to reproduce experimental (ro)vibrational energies to within tenths of cm(-1) if they are included in the refinement or interpolated between states that are involved in the optimization, whereas extrapolated states tend to have somewhat larger remaining discrepancies; (4) the PES refined only for the j = 0 states introduces a minor systematic error for J> 0 states, with discrepancies growing with j; (5) when the number of experimental energies included in the refinement greatly exceeds the number of grid points upon which the PES is optimized, the systematic error of treating states with different J rotational quantum numbers can be reduced and an impressive average accuracy can be achieved for all rovibrational states; and (6) in the case of quasibound (also known as resonance) rovibrational states, energies can be computed to accuracies similar to those of the bound states and excellent lifetimes (widths) can also be determined. Changes in thermochemical functions upon inclusion of quasibound states during direct summation is discussed.