Nonlinear gyrokinetic simulations are conducted to investigate turbulent transport in tokamak plasmas with rotational shear. At sufficiently large flow shears, linear instabilities are suppressed, but transiently growing modes drive subcritical turbulence whose amplitude increases with flow shear. This leads to a local minimum in the heat flux, indicating an optimal E x B shear value for plasma confinement. Local maxima in the momentum fluxes are observed, implying the possibility of bifurcations in the E x B shear. The critical temperature gradient for the onset of turbulence increases with flow shear at low flow shears; at higher flow shears, the dependence of heat flux on temperature gradient becomes less stiff. The turbulent Prandtl number is found to be largely independent of temperature and flow gradients, with a value close to unity.