$\newcommand\F{\mathbf{F}}$
$\newcommand\Z{\mathbf{Z}}$
$\newcommand\SL{\mathrm{SL}}$
Because of the existence of the Weil paring,
elliptic curves with such a subgroup only exist when
$p \equiv 1 \mod \ N$.
Let $S_N$ denote the set of elliptic curves over $\F_p$ such
that $E[N]$ is defined over $\F_p$.
It will be slightly easier to assume that $N \ge 3$. In this case,
$Y(N)$ is a fine moduli space, and an $\F_p$-point on $Y(N)$
corresponds to
a pair $(E,\alpha:E[N] \simeq \Z/N\Z \times \Z/N \Z)$ defined over
$\F_p$. Given an elliptic curve $E \in S_N$,
how many points does it contribute to $Y(N)$? For a curve $E$
whose automorphism group is $\Z/2\Z$, We see that out of
the $|\SL_2(\Z/N\Z)|$ possible choices of $\alpha$
(technical remark, we have fixed a Weil pairing so that $Y(N)$ is
connected), $(E,\alpha) \simeq (E,\alpha')$ only if $\alpha' = \alpha$
or $\alpha' = [-1] \alpha$. Thus
$E$ contributes $|\SL_2(\Z/N\Z)|/2$ points to $Y(N)(\F_p)$.
In general, $E$ may have slightly more
automorphisms, and we deduce that (for $N \ge 3$):
$$|Y(N)(\F_p)| = |\SL_2(\Z/N\Z)| \sum_{E \in S_N}
\frac{1}{|\mathrm{Aut}(E)|}.$$
Note that the quantity on the right is very close to
$|\SL_2(\Z/N\Z)| \cdot |S_N|/2$, one only has to worry about the
elliptic curves with $j = 0$ or $j = 1728$, and this can be done by hand
if one wants to cross all the i's and dot all the t's.
Suppose that $X(N)$ has $c_N$ cusps and genus $g_N$ (there are some
explicit slightly unpleasant formulas for these numbers, which can
be found (for example) in Shimura's book. All the cusps
are defined over $\F_p$ (with $p \equiv 1 \mod N$) so by the
Riemann hypothesis for finite fields,
$$|Y(N)(\F_p) - (1+p) + c_N| = |X(N)(\F_p) - (1+p)| \le
2 g_N \sqrt{p}.$$
If $g_N = 0$ (which only happens if $N \le 5$), this leads to an
exact formula for $|S_N|$.
In general, at least for large $p$, we see that
$$|S_N| \sim \frac{2p}{|\SL_2(\Z/N \Z)|}.$$
To make this all completely explicit for $N = 3$ (for example),
one gets, presuming I have not made a horrible computational error which is quite possible:
$$S_3 = \begin{cases} (p+11)/12, & p \equiv 1 \mod \ 12, \\\
(p+5)/12, & p \equiv 7 \mod \ 12. \end{cases}$$
(note that $p \equiv 1 \mod 3$):
Of course, "exact formulas" will only exist for $N \le 5$.
Some related and slightly more difficult counting problems are also
nicely explained by Lenstra here (See 1.10):
https://openaccess.leidenuniv.nl/bitstream/1887/3826/1/346_086.pdf
Best Answer
A standard method to find such an equation (in principle) is the following, shown to me by Tate, but maybe due to Mordell(?). If $P_0$ is a point of order at least $4$ on $E$ (infinite order is allowed), then after a change of coordinates, we can find an equation for $E$ of the form $$ E : y^2 + u x y + v y = x^3 + v x^2 \quad\text{with}\quad P_0=(0,0). $$ (It's a nice exercise to check this assertion.) So now we start with this equation and point, and we want $P_0$ to have order $N$. If $N$ is odd, say $N=2n+1$, then use the group law on $E$ to compute $x(nP_0)$ and $x\bigl((n+1)P_0\bigr)$. These will each be rational functions of $u$ and $v$; i.e., they're in $\mathbb Q(u,v)$. Setting $$ x(nP_0)=x\bigl((n+1)P_0\bigr)$$ and clearing the denominators will give you a polynomial in $(u,v)$ whose vanishing is an affine plane model for $X_1(N)$ if $N$ is prime. However it tends to be rather singular at points where the discriminant of $E$ vanishes. And if $N$ is composite, then the polynomial will factor, i.e., the curve will have components for each divisor of $N$ greater than $3$. Finally, if $N=2n$ is even, you can do the same thing with the decomposition $N=(n+1)+(n-1)$, or you could set the numerator of $(2y+ux+v)(nP_0)$ equal to $0$, which is another way of forcing $2nP_0$ to be $O$.