Results from a fully nonlinear three-dimensional toroidal electrostatic gyrokinetic simulation of the ion temperature gradient instability are presented. The model has adiabatic electrons and the complete gyrophase-averaged ion dynamics, including trapped particles. Results include the confirmation of the radially elongated ballooning mode structure predicted by linear theory, and the nonlinear saturation of these toroidal modes. The ensuing turbulent spectrum retains remnants of the linear mode structure, and has very similar features as recent experimental fluctuation measurements.PACS numbers: 52.35. Ra, 52.25.Gj, 52.35.Qz Recent advances in both nonlinear df methods for gyrokinetic simulation [1,2] and massively parallel supercomputing now make it possible to simulate a sizable fraction of a tokamak plasma using realistic physical parameters. Here, we report results from the first whole cross section three-dimensional (3D) electrostatic toroidal gyrokinetic simulation. We investigate the nonlinear evolution of the ion temperature gradient (ITG) driven instability and the associated turbulence and transport in realistic geometry and dimensionality. The ITG mode has long been considered a plausible candidate to explain the observed anomalous ion heat transport in tokamak plasmas, which is substantially above the predicted neoclassical values [3,4]. The simulation results presented below show very similar features in terms of the fluctuation spectrum as the recent beam emission spectroscopy (BES) diagnostic on TFTR [5].In these simulations, the ions are fully gyrokinetic, including trapped particles. The electrons are treated as adiabatic which permits a moderate size time step (simulations with kinetic electrons are feasible, but the time step would need to be smaller by the factor ivAv/). The simulation is running efficiently on massively parallel supercomputers (currently the CM-200 and CM-5) which allow simulations of relatively large systems (e.g., a ^ lOOp/ minor radius. Ax ^=^ pi). Typical runs up to this point have used 10^ to 10^ particles with 1 to 2 particles per grid cell, and a timing of 2-3 jUS per particle per time step on a fully configured CM-200. Fine grid resolution is needed in the toroidal direction because the mode structure is helical (elongated along the magnetic field lines, i.e., k\\<^k±), resulting in a smaller number of particles per grid cell relative to conventional slab simulations.Starting with the electrostatic gyrokinetic equations with a nonuniform equilibrium B field [6], we write /(z,r) ==/o(z) + ^/(z,/), where z = (R,rii,/x), and expand z into its equilibrium and perturbed parts: z=z^-f-z^ /o(z) is a Maxwellian and satisfies z^S^o(z)=0. The equation for the perturbed ion distribution function 5f is then [1] S,5/+z9z^/=--z^9j^o,where the magnetic moment fi is time independent and the other equilibrium and perturbed phase space variables are evolved using (R",r,r)=(-riib*-h^bxV5,b*-A/V5), tLu^ ii /,i^ = (R\n|)