A new gyrokinetic toroidal particle code has been developed to study the ion temperature gradient (ITG) driven turbulence in reactor relevant tokamak parameters. We use a new method based on a canonical Maxwellian distribution F-CM(Pphi, epsilon, mu), which is defined by three constants of motion in the axisymmetric toroidal system-the canonical angular momentum Pphi, the energy epsilon, and the magnetic moment mu. A quasi-ballooning representation enables linear and nonlinear high-m, n global calculations to be carried out, with a good numerical convergence. Conservation properties are improved by using optimized particle loading. From comprehensive linear global analyses over a wide range of unstable toroidal mode numbers (n = 0-100) in large tokamak parameters (a/rho(ti) = 320-460), it is found that the reversed shear configuration produces an effective stabilizing effect on the ITG mode in the q(min) region through global effects. In the nonlinear simulation, it is found that the new method based on F-CM can simulate a zonal flow damping correctly; and spurious zonal flow oscillations, which are observed in a conventional method based on a local Maxwellian distribution F-LM(psi, epsilon, mu), do not appear in the nonlinear regime.