  We present a new formulation for numerically obtaining axisymmetric
equilibrium structures of rotating stars in two spatial dimensions. With a view
to apply it to the secular evolution of rotating stars, we base it on the
Lagrangian description, i.e., we solve the force-balance equations to find the
spatial positions of fluid elements endowed individually with a mass, specific
entropy and angular momentum. The system of nonlinear equations obtained by
finite-differencing the basic equations are solved with the W4 method, which is
a new multi-dimensional root-finding scheme of our own devising. We augment it
with a remapping scheme to avoid distortions of the Lagrangian coordinates. In
this first one of a series of papers, we will give a detailed description of
these methods initially. We then present the results of some test calculations,
which include the construction of both rapidly rotating barotropic and
baroclinic equilibrium states. We gauge their accuracies quantitatively with
some diagnostic quantities as well as via comparisons with the counterparts
obtained with an Eulerian code. For a demonstrative purpose, we apply the code
to a toy-model cooling calculation of a rotating white dwarf.
