  Multiple zeta values (MZVs) are generalizations of Riemann zeta values at
positive integers to multiple variable setting. These values can be further
generalized to level $N$ multiple polylog values by evaluating multiple
polylogs at $N$-th roots of unity. In this paper, we consider another level $N$
generalization by restricting the indices in the iterated sums defining MZVs to
congruences classes modulo $N$, which we call the MZVs at level $N$. The goals
of this paper are two-fold. First, we shall lay down the theoretical
foundations of these values such as their regularizations and double shuffle
relations. Second, we will generalize the multiple divisor functions (MDFs)
defined by Bachman and K\"uhn to arbitrary level $N$ and study their relations
to MZVs at level $N$. These functions are all $q$-series and similar to MZVs,
they have both weight and depth filtrations. But unlike that of MZVs, the
product of MDFs usually has mixed weights; however, after projecting to the
highest weight we can obtain an algebra homomorphism from MDFs to MZVs.
Moreover, the image of the derivation $\mathfrak{D}=q\frac{d}{dq}$ on MDFs
vanishes on the MZV side, which gives rise to many nontrivial
$\mathbb{Q}$-linear relations. In a sequel to this paper, we plan to
investigate the nature of these relations.
