forms for $\mu$. All previous works (@Tomimura_Eriguchi_2005, @Yoshida_Eriguchi_2006, @Yoshida_Yoshida_Eriguchi_2006, @Lander_Jones_2009, @Otani_Takahashi_Eriguchi_2009) fixed $\mu = \mu_0 $ (constant). We choose a different functional form for $\mu$ in this paper as follows: $$\begin{aligned}
   \mu(\Psi)                & = & \mu_0 (\Psi + \epsilon )^m \ , \\
   \int \mu(\Psi) \, d \Psi & = & \frac{\mu_0}{m+1}(\Psi + \epsilon )^{m+1} \ ,
   \label{eq:muuu}\end{aligned}$$ where $m$ and $\epsilon$ are two constant parameters. In order to avoid singular behavior, we fix $\epsilon = 1.0 \times 10^{-6}$ in all calculations. As we shall see below, the parameter $m$ determines a degree of localization of the interior [*poloidal*]{} magnetic field. We assume that [*poloidal*]{} magnetic fields extend throughout the whole space and that there are no discontinuities even at the stellar surface. The global magnetic field configurations of our models are nearly dipole-like because of the requirement of the functional form for $\kappa$ at the stellar surface. These configurations contain closed [*poloidal*]{} magnetic field lines inside the star. The flux function $\Psi$ attains its maximum at the central parts of these closed field lines and it takes its minimum on the symmetric axis and at infinity. The minimum value is zero because of $\Psi = r \sin \theta A_\varphi$ 