  We develop a theory connecting the following three areas: (a) the mean field
equation (MFE) $\triangle u + e^u = \rho\, \delta_0$, $\rho \in \mathbb R_{>0}$
on flat tori $E_\tau = \mathbb C/(\mathbb Z + \mathbb Z\tau)$, (b) the
classical Lam\'e equations and (c) modular forms. A major theme in part I is a
classification of developing maps $f$ attached to solutions $u$ of the mean
field equation according to the type of transformation laws (or monodromy) with
respect to $\Lambda$ satisfied by $f$.
  We are especially interested in the case when the parameter $\rho$ in the
mean field equation is an integer multiple of $4\pi$. In the case when $\rho =
4\pi(2n + 1)$ for a non-negative integer $n$, we prove that the number of
solutions is $n + 1$ except for a finite number of conformal isomorphism
classes of flat tori, and we give a family of polynomials which characterizes
the developing maps for solutions of mean field equations through the
configuration of their zeros and poles. Modular forms appear naturally already
in the simplest situation when $\,\rho=4\pi$.
  In the case when $\rho = 8\pi n$ for a positive integer $n$, the solvability
of the MFE depends on the \emph{moduli} of the flat tori $E_\tau$ and leads
naturally to a hyperelliptic curve $\bar X_n=\bar X_{n}(\tau)$ arising from the
Hermite-Halphen ansatz solutions of Lam\'e's differential equation $\frac{d^2
w}{dz^2}-(n(n+1)\wp(z;\Lambda_{\tau}) + B) w=0$. We analyse the curve $\bar
X_n$ from both the analytic and the algebraic perspective, including its local
coordinate near the point at infinity, which turns out to be a smooth point of
$\bar{X}_n$. We also specify the role of the branch points of the hyperelliptic
projection $\bar X_n \to \mathbb P^1$ when the parameter $\rho$ varies in a
neighborhood of $\rho = 8\pi n$.
