It is now widely believed that low frequency turbulence developing from smallscale instabilities is responsible for the phenomenon of anomalous transport generally observed in magnetic confinement fusion experiments. The microinstabilities are driven by gradients of equilibrium density, ion and electron temperatures and magnetic field strength. Gyrokinetic theory is based on the Vlasov-Maxwell equations and, consistent with the ordering, averages out the fast particle gyromotion, reducing the phase space from 6 to 5 dimensions. Solving the resulting equations is a non-trivial task. Difficulties are associated with the magnetic confinement geometry, the strong disparities in space and time scales perpendicular and parallel to B, the different time scales of ion and electron dynamics, and the complex nonlinear behaviour of the system. The main numerical methods are briefly presented together with some recent developments and improvements to the basic algorithms. Recent results are shown, with emphasis on the roles of zonal E x B flows, of parallel nonlinearity and of toroidal coupling on the saturation of ion temperature gradient (ITG) driven turbulence in tokamaks.