We derive a system of equations to describe stationary axisymmetric MHD equilibria characterized by toroidal and poloidal flows as well as plasma anisotropy due to strong magnetic field and eventual auxiliary heating methods. The system consists of a nonlinear partial differential equation for the poloidal magnetic flux function and an algebraic Bernoulli-type equation, which relates plasma density with several surface functions. We analyse the ellipticity of the equation and the plasma density bifurcation as the Alfvén-Mach number is changed.