Understanding the turbulent dynamics of the tokamak scrape-off layer (SOL) is critical for the safe operation of future fusion reactors. In this region of the tokamak device, the plasma exhausts particles and heat onto plasma facing hardware components, which may limit their lifetime if material constraints are not met. Recently, we have gained deep insights into SOL dynamics by means of electromagnetic fluid turbulence simulations carried out with the Global Braginskii Solver (GBS). GBS is currently capable of carrying out simulations of the SOL of medium-size tokamaks such as TCV or Alcator C-Mod. We present the new version of the GBS code, which is intended to address larger SOL plasmas using advanced two-fluid models. GBS has undergone major transformations, including (a) porting of the code kernel to F2008, (b) upgrading pure MPI parallelism to employ a 3D cartesian communicator, and (c) including a matrix-free parallel multigrid solver to address generalized Poisson and Ampere operators. There is also an ongoing effort to port this code to MIC/GPU architecture. The new GBS code achieves excellent parallel scalability up to 8192 cores for a C-Mod size plasma (e.g. weak/strong scalings), while at the same time it is possible to achieve a fast and scalable solution of time dependent operators. As an example, we discuss two new physics capabilities of GBS allowed by numerical advances. The first is the inclusion of electromagnetic fluctuations at realistic plasma sizes, which is now possible due to the reduced numerical cost and improved parallel scalability. Second, the Boussinesq approximation in the evaluation of the divergence of the polarization drift is not justified in the tokamak scrape-off layer, where dn/n ~ 1. The drift-reduced equations have been rederived, and we have implemented and benchmarked a non-Boussinesq Poisson operator in GBS. Finally, we present the latest GBS simulation results using realistic size and parameters from inner-wall limited discharges, which is the chosen scenario for ITER start-up. Our simulations reveal the presence of steep plasma gradients near the last closed flux surface, which implies that a large portion of the exhaust heat can be deposited in a very small area. The turbulent dynamics recovered show low frequency, large amplitude drift and ballooning modes with a 1cm radial and poloidal correlation lengths. Intermittent transport events are present, leading to strongly skewed fluctuation probability distribution functions. The simulations have been benchmarked against measurements from the same experiments, showing very good agreement in many observables.