Let me try to explain why this problem is hard (concurring with Alexey Ustinov here). Fourier analysis shows that $$ N_p = \frac{1}{p^2} \sum_{m=0}^{p-1} \sum_{n=0}^{p-1} \left( \sum_{x=1}^{p-1} e_p( nx + mx^{-1} )\right) \left(\sum_{x=1}^{(p-1)/2} e_p(-nx) \right) \left(\sum_{x=1}^{(p-1)/2} e_p(-mx) \right) $$
where $ \left( \sum_{x=1}^{p-1} e_p( nx + mx^{-1} \right)$ is a Kloosterman sum $K(mn;c)$ and $ \left(\sum_{x=1}^{(p-1)/2} e_p(-nx) \right) \left(\sum_{x=1}^{(p-1)/2} e_p(-mx) \right)$ is large for $n$, $m$ close to zero and small elsewhere.
The term $n=0, m=0$ gives the main term, where $n=0$ and $m \neq 0$ or vice versa gives a small and easily computable term, and the interesting error term comes from the range where $n$ and $m$ are nonzero. This can be represented as a linear combination of Kloosterman sums $K(mn;p)$, with the largest contributions coming from $mn = \pm 1$.
We expect the sums $K(a;p)/\sqrt{p} $ for any $a$ to be equidistributed as random variables with a Sato-Tate (semicircular) distribution. It's reasonable to also conjecture independence over different $a$, which would enable calculating the distribution of any finite sum of Kloosterman sums. (In fact this follows from a conjecture I made in my PhD thesis). This is an infinite sum of Kloosterman sums, but it is probably reasonable to make a similar independence statement, which should give a Gaussian distribution or one very similar to it.
However, almost nothing is known about the distribution of $K(a;p)$ for fixed $a$ and varying $p$. (The case of fixed $p$ and varying $a$ is well-known by work of Deligne and Katz). Even the first moment is not known, although we do have a bound for the average of $K(a;c)$ over all numbers $c$, proved using the Kusnetsov formula.
This sum of Kloosterman sums looks much harder to analyze than the distribution of a single Kloosterman sum, so it is probably also too hard for current techniques.