An efficient iterative, free-plasma-boundary solver for the Grad–Shafranov–Bernoulli system of equations, that describes the ideal MHD equilibrium of a toroidally axisymmetric plasma with flow, is presented. The code implements a numerical scheme recently developed in the context of free-plasma-boundary solvers for ideal static MHD equilibria with magnetic islands and stochastic regions for stellarators. The shape of the plasma edge is permitted to change as needed until the total net force eventually vanishes en route to the equilibrium. Complex coil configurations can be treated in the toroidally axisymmetric approximation. The code opens the possibility of quantifying the changes that plasma flows may induce on important features of a tokamak equilibrium such as the shape of the plasma edge, the plasma confining volume, the position of the magnetic axis or the position of the X-point, among others. Some examples, selected for illustrative purposes, are shown for the ITER baseline magnetic configuration.