We present a framework for the preconditioned iterative solution of systems arising from variational mean field games, with local couplings and either periodic or Neumann boundary conditions, that allows the application of different parallelization strategies across the time variable. Following a variational approach, we employ a finite difference discretization and solve the resulting finite-dimensional optimization problem using the Chambolle--Pock primal--dual algorithm. As this requires the repeated solution of linear systems within the evaluation of proximal operators, we embed within our solver a general class of parallel-in-time preconditioners based on diagonalization techniques in the temporal direction, implemented via discrete Fourier transforms. These enable efficient, scalable iterative solvers for each linear system, with robustness across a wide range of viscosities. For structured grids, we further develop fast recursive solvers that exploit tensor-product structure, while retaining flexibility for more general geometries. The key contributions of this paper includes the practical and robust numerical solution of very large mean field games systems, previously untested in the literature to our knowledge, in particular systems within the numerically challenging regime of small but non-zero viscosities. Numerical experiments demonstrate that the proposed preconditioning strategies significantly reduce iteration counts and wall-clock time, and exhibit favorable parallel scalability.