The feasibility of using full ion kinetics, instead of gyrokinetics, in simulating low-frequency Ion-Temperature-Gradient (ITG) instabilities in tokamaks has recently been demonstrated. The present work extends the full ion kinetics to the nonlinear regime and investigates the nonlinear saturation of a single-n ITG instability due to the E × B trapping mechanism (n is the toroidal mode number). The saturation amplitude predicted by the E × B trapping theory is found to agree with the saturation level observed in the simulation. In extending to the nonlinear regime, we developed a toroidal Boris full orbit integrator, which proved to be accurate in capturing both the short-time scale cyclotron motion and long time scale drift motion, with good kinetic energy conservation and toroidal angular momentum conservation in tokamak equilibrium magnetic fields. This work also extends the previous work from analytic circular magnetic equilibria to general numerical magnetic equilibria, enabling simulation of realistic equilibria reconstructed from tokamak experiments.