The acceleration at the point of reflection is actually quite complicated. It is caused by the elastic forces of the surface and the ball and has a complicated time dependence. However, the timespan in which the ball touches the ground is very short (especially if the ground and the ball are very rigid), therefore we can simplify the actual acceleration curve for kinematic calculations and simply assume that the velocity jumps instantaneously from $v$ to $-v$ in the moment of contact (because the contact with the ground during which this acceleration acts is perhaps on the order of a few millisecons, while we are interested in time scales $\Delta t$ and corresponding length scales $v \Delta t$ a range of seconds).
In this approximation, which is identical to taking the limit of infinite rigidity of ball and surface (formally described by the bulk modulus $E$), the acceleration of the ball is given by
$$ a(t) = -g\vec e_z - 2 v(t) \delta(t - t^*), $$
where $t^*$ is the moment of contact.
Here $\delta(t)$ is the so called Dirac delta function, which can be though of as a very narrow but high peak whose integral is one. There are several ways to define the delta function, the most useful one for physicists is the rule
$$ \int_{-\infty}^\infty dx\, \delta(x - y) f(y) = f(x) $$
for all functions $f$ (if you do this properly the functions must have some regularity). You can also always think of the $\delta$-function as the limit of narrow peak functions (where you can use any family of peak function you like that have unit area and concentrate more and more around zero):
$$ \delta(x) = \lim_{\eta \to 0} \frac{1}{\sqrt{2\pi \eta}} e^{-x^2/\eta} $$
to calculate something rigorously, you have to take the limit after integration (then it will often exist formally, while the $\delta$-function does not exists formally as a limiting function).
Another useful representation of the $\delta$-function comes from the theory of Fourer transformation:
$$ \delta(x) = \frac{1}{2\pi} \int_{-\infty}^\infty dk\, e^{ikx}. $$
With those properties we can simply check that the above acceleration profile gives the result for $v(t)$:
\begin{align*}
v(t) &= \int_0^t dt\, a(t) = - g t - \int_0^t dt\, 2 v(t) \delta(t - t^*) \\
&= - g t + 2 g t^* \theta(t - t^*).
\end{align*}
Where
$$\theta(t) = \begin{cases} 0 & t < 0 \\ 1 & t > 0 \end{cases} $$
is the Heavyside step function (in short $\theta$-function).
$z(t)$ is given by:
$$z(t) = z(0) + \int_0^t v(t) = -\frac 1 2 gt^2 + \theta(t - t^*) 2 g t^* t. $$
Of course, one can proceed to express $t^*$ in terms of the initial height.
