A four-field system of equations which includes the neoclassical flow damping effects and the lowest-order finite-Larmor-radius (FLR) corrections is deduced from a system of gyrofluid equations. The FLR corrections to the poloidal flow damping are calculated by solving a simplified version of the gyrokinetic equation. This system of equations is applied to the study of a chain of freely rotating magnetic islands in a tokamak, resulting from the nonlinear evolution of a resistive tearing mode, to determine the island rotation velocity consistently with the fields' radial profiles close to the resonant surface. The island rotation velocity is determined by imposing the torque balance condition. The equations thus deduced are applied to the study of two different collisionality regimes, namely the weak-damping regime and the intermediate-damping regime. The equations reduce, in the weak-damping regime, to a form already obtained in previous works, while an additional term, containing the lowest order FLR corrections to the poloidal flow damping, appears in the intermediate-damping regime. The numerical integration of the final system of equations allows the determination of the dependence of the island rotation velocity on the plasma collisionality and the island width compared to the ion Larmor radius. The results show that, in the intermediate-damping regime, the island rotation velocity is almost completely determined by the neoclassical effects, with the island width playing a minor role. The parameter ηi=Ln/LT, where Ln and LT are the density and temperature gradient length scales, plays an important role in determining the island rotation velocity.