Here is a solution for $\frac{p-1}{2}$ of the level sets. Let $\varphi=x^2+y^2-x+y-axy$, and assume that $a\not=\pm 2$.
Following the suggestion in https://mathoverflow.net/q/356936 we first apply the change of variables given by $x= z+(a/2)y$ to obtain $$\varphi=z^2-z-b(y^2+(1+a/2)^{-1}y),$$ where $b=a^2/4-1$. Then letting $u=z-1/2$, $v=y+(2+a)^{-1}$, and $c=\frac{-1}{4}\left(\frac{1-a/2}{1+a/2}\right)$ we have $$\varphi=u^2-bv^2+c.$$
To compute the level sets you need to compute the fibres $\varphi^{-1}(d)$ for an image point $d$. If $d-c$ is a quadratic residue, then $d-c=e^2$ and so letting $s=u/e$ and $t=v/e$ we obtain: $$s^2-bt^2=1.$$
Having all solutions $(s,t)$ determines the solutions $(u,v)$, in turn the solutions $(y,z)$, and in turn the solutions $(x,y)$.
So it suffices to find the solutions to $s^2-bt^2=1.$ But such solutions are explicitly determined in:
Tekcan, Ahmet, The number of rational points on conics $C_{p,k}:x^2βky^2=1$ over finite fields $\mathbb{F}_p$. Int. J. Math. Sci. 1 (2007), no. 2, 150β153.
MathSciNet Summary:
"Let $p$ be a prime number, $\mathbb{F}_p$ be a finite field, and let
$k\in \mathbb{F}_p^*$. In this paper, we consider the number of
rational points on conics $C_{p,k}: x^2βky^2=1$ over $\mathbb{F}_p$.
We prove that the order of $C_{p,k}$ over $\mathbb{F}_p$ is $pβ1$ if
$k$ is a quadratic residue mod $p$ and is $p+1$ if $k$ is not a
quadratic residue mod $p$....''
I am not sure what happens when $d-c$ is a non-residue (I haven't tried), but since $\varphi$ is onto and $(p-1)/2$ of the non-zero values are quadratic residues, this gets you roughly half way.