Intrinsic toroidal plasma rotation due to the neoclassical toroidal plasma viscosity (NTV) effect induced by a three-dimensional helical magnetic field ripple in tokamaks is investigated in this paper. The intrinsic rotation is determined self-consistently by searching for the roots of the ambipolarity constraint, after evaluation of the particle fluxes from the numerical modelling. In the low-collisionality case, there are three roots, in which two are stable roots. One corresponds to the 'ion root' in the counter-current direction, and the other stable one corresponds to the 'electron root' in the co-current direction, near which the electron flux is dominant. Both of the two stable roots scale like the diamagnetic frequency. In the high-collisionality case, there is only one 'ion' root. The application of this modelling for International Thermonuclear Experimental Reactor (ITER) cases is discussed. In a large range of plasma radii, there are three roots. The NTV torque drives plasma rotation in ITER towards one of the stable roots, depending on the initial condition. The amplitudes of the electron roots near the pedestal in both baseline and steady-state scenarios are much larger than that of the ion roots. The amplitudes of the NTV torque density and the electron roots near the pedestal increase with increasing height of the temperature pedestal in the ITER baseline scenario.