For example, suppose you just have the one equation, say $x^2 + 1 = 0$. The
Grobner basis will consist of $x^2+1$. This has no bearing on whether $-1$ is a square in your field.

EDIT: For finding solutions in your field, the following **might** work (in particular if the ideal is zero-dimensional).

Take a "plex" Gröbner basis, and factor the first polynomial in it over your field. Note that factoring a polynomial is pretty efficient, even modulo a $128$-bit prime. Any factors that are linear in one of the variables (say $x$) correspond to candidates for values of $x$, possibly dependent on other variables. For each of these,
substitute that value for $x$ into the rest of your basis and continue. Any factors that are nonlinear and contain just one variable can be disregarded, as they won't correspond to solutions in the field. Factors with more than one variable, and no variable for which the factor has degree $1$, might be problematic, though.

For example, consider the polynomials
$$ \matrix{21\,xy-121\,{z}^{2}+6\,x+7\,y-110\,z-23 \cr 33\,xz-49\,{y}^{2}+18\,x-42\,
y+22\,z+3\cr-9\,{x}^{2}+77\,zy-6\,x+42\,y+11\,z+6 \cr {z}^{4}+{y}^{3}+{x}^{2
}+230208773171659848188232517051917004253}
$$
over $F_p$ where $p$ is the prime $850705917302346158658436518579420528641$.

In Maple:

```
> with(PolynomialIdeals):
> p:= 850705917302346158658436518579420528641;
> G := [21*x*y-121*z^2+6*x+7*y-110*z-23, 33*x*z-49*y^2+18*x-42*y+22*z+3,
> -9*x^2+77*y*z-6*x+42*y+11*z+6,
> z^4+y^3+x^2+230208773171659848188232517051917004253];
> J:= PolynomialIdeal(G, characteristic=p);
> B:= Groebner[Basis](J, plex(x,y,z));
> B:= map(Factor, B) mod p;
```

$$B:= [ \left( {z}^{3}+638987694045812089487943040989684012240\,{z}^{2}+
788559908137201150047779839724942410702\,z+
166287164820212048829798599861338659586 \right) \left( z+
154673803145881119756079367014440096117 \right) ,y+
121529416757478022665490931225631504090\,z+
729176500544868135992945587353789024549,
567137278201564105772291012386280352426+
283568639100782052886145506193140176210\,z+x]
$$

The first polynomial depends only on $z$ and there is one linear factor.
Substituting the corresponding value $z = 696032114156465038902357151564980432524$ into the other two polynomials in the Gröbner basis, you get corresponding values for $x$ and $y$. There is just this one solution of the equations over $F_p$.