  (Abridged) We develop a theoretical framework to construct axisymmetric
magnetic equilibria in stars, consisting of both poloidal and toroidal magnetic
field components. In a stationary axisymmetric configuration, the poloidal
current is a function of the poloidal magnetic flux only, and thus should
vanish on field lines extending outside of the star. Non-zero poloidal current
is limited to a set of toroid-shape flux surfaces fully enclosed inside the
star. If we demand that there are no current sheets, then on the separatrix
delineating the regions of zero and finite toroidal magnetic field both the
poloidal flux function and its derivative should match. Thus, for a given
magnetic field in the bulk of the star, the elliptical Grad-Shafranov equation
that describes magnetic field structure inside the toroid is an ill-posed
problem, with both Dirichlet and Newman boundary conditions and {\it a priori}
unknown distribution of toroidal and poloidal electric currents. We discuss a
procedure which allows to solve this ill-posed problem by adjusting the unknown
current functions. We find a poloidal current-carrying solution that leaves the
shape of the flux function and, correspondingly, the toroidal component of the
electric current the same as in the case of no poloidal current. The equilibria
discussed in this paper may have arbitrary large toroidal magnetic field, and
may include a set of stable equilibria. The method developed here can also be
applied to magnetic structure of differentially rotating stars, as well as to
calculate velocity field in incompressible isolated fluid vortex with a swirl.
