The melting-like transition in sodium clusters Na-N, with N = 55, 92, and 132, is studied by using constant-energy molecular dynamics simulations. An orbital-free version of the Car-Parrinello technique is used that scales linearly with system size, allowing for investigations of the thermal behavior of large clusters. The ground-state isomer of Na-142 tan incomplete three-shell icosahedron) melts in two steps: the first (at similar to 240 K) is characterized by the high mobility of atoms located on the cluster surface, and the second, homogeneous melting (at similar to 270 K) involves diffusive motion of all of the atoms across the cluster. For the case of Na-92, the icosahedral structure has a larger number of surface vacancies and melts in two well-separated steps, with surface melting at similar to 130 K and homogeneous melting at similar to 240 K. Na-55, a complete two-shell icosahedron, melts in a single stage at similar to 190 K. Our results on homogeneous melting for Na-142 and Na-92 are in excellent agreement with recent experimental determinations of melting temperatures and latent heats. However, the experimentally observed enhancement of the melting temperature around N = 55 is not reproduced by the calculations.