We shall use some ideas and observations from the answer at Controlling subsolutions of a second order linear ODE, as well as some additional considerations to provide an answer to the above question.
In fact, we shall provide the exact upper bound on $y$ -- for each given quadruple of values $(a,b,f_0,g_0)$.
It was shown in the just mentioned answer that, for any initial conditions on $f(0)$ and $f'(0)$, we have $f(x_n)\le0$ for some sequence $(x_n)$ converging to $\infty$ and all natural $n$.
Therefore, there indeed exists the real number $y$ in question. Moreover, we have $y\in(0,\infty)$, $f(y)=0$, and $f>0$ on the interval $(0,y)$.
As in the mentioned answer, introduce
\begin{equation}
r(x):=f(x)/e^{ax}. \tag{1}
\end{equation}
Then inequality (*) in the above question can be rewritten as
\begin{equation}
r''+b^2r\le0, \tag{2}
\end{equation}
and we also have $r>0$ on the interval $(0,y)$, with $r(y)=0$.
So, $r''\le-b^2r<0$ and hence $r$ is strictly concave on $(0,y)$.
If $f(0) = f_0\in(0,\infty)$, then $r(0) = f_0\in(0,\infty)$. Alternatively, if $f(0) = 0$ and $f'(0) = g_0\in(0,\infty)$, then $r(0) = f_0=0$ and $r'(0)[= g_0-af_0] = g_0\in(0,\infty)$.
So, one of the following two cases must take place:
Case 1: $r(0) = r_0:=f_0\in[0,\infty)$ and $r'(0)=g_{a,0}:=g_0-af_0\in(0,\infty)$, or
Case 2: $r(0) = r_0\in(0,\infty)$ and $r'(0)=g_{a,0}\in(-\infty,0]$.
Consider now Case 1. Since $r''<0$ on $(0,y)$ and $r'(0)>0$, we have $r'(y)<0$ -- because otherwise we would have $r'>0$ on $(0,y)$ and hence $0=r(y)>r(0)\ge0$, a contradiction. So,
there exists a unique $t\in(0,y)$ such that $r'(x)$ (strictly) decreases from $r'(0)>0$ to $r'(t)=0$ to $r'(y)<0$ as $x$ increases from $0$ to $t$ to $y$. So, $r'>0$ on $(0,t)$ and $r'<0$ on $(t,y)$, which in turn implies that $r$ continuously increases from $r(0)=r_0\in[0,\infty)$ to
\begin{equation*}
r_{\max}:=r(t)>r_0\ge0
\end{equation*}
on $[0,t]$, and $r$ continuously decreases from $r_{\max}>0$ to $r(y)=0$ on $[t,y]$.
So, for the restrictions $\rho_1:=r|_{I_1}$ and $\rho_2:=r|_{I_2}$ of the function $r$ to each of the intervals $I_1:=[0,t]$ and $I_2:=[t,y]$, the corresponding inverse functions $\rho_j^{-1}$ are well defined, and then so are the functions $p_j:=r'\circ\rho_j^{-1}$,
so that $r'(x)=p_j(\rho_j(x))=p_j(r(x))$ for $x\in I_j$, for each $j=1,2$.
By the chain rule, $r''=p_j'(r)r'=p_j'(r)p_j(r)=q_j'(r)/2$, where
$$q_j(r):=p_j(r)^2,$$
twice the "kinetic energy".
Thus, inequality (2) can be rewritten as
\begin{equation}
q_j'(r)\le-2b^2r, \tag{3}
\end{equation}
for $r$ in the corresponding the intervals $R_1:=[r_0,r_{\max}]$ and $R_2:=[0,r_{\max}]$, which in turn correspond to the intervals $I_1=[0,t]$ and $I_2=[t,y]$ for the values of $x$.
Integrating (3), we have
\begin{equation}
{r'}^2=p_1(r)^2=q_1(r)\le q_1(r_0)-b^2(r^2-r_0^2)=g_{a,0}^2-b^2(r^2-r_0^2) \tag{4}
\end{equation}
for $r\in[r_0,r_{\max}]$, that is, for $x\in[0,t]$.
Since $r'(t)=0$ and $r(t)=r_{\max}$, it follows from (4) that
\begin{equation}
r_{\max}\le r_*:=\sqrt{r_0^2+g_{a,0}^2/b^2}. \tag{5}
\end{equation}
Another integration of (3) yields
\begin{equation}
-{r'}^2=-q_j(r)=q_j(r_{\max})-q_j(r)\le-b^2(r_{\max}^2-r^2) \tag{6}
\end{equation}
for $r\in R_j$, that is, for $x\in I_j$.
Since $r'\ge0$ on $I_1=[0,t]$ and $r'\le0$ on $I_2=[t,y]$, (6) can be rewritten as
\begin{equation}
\text{$r'\ge b\sqrt{r_{\max}^2-r^2}$ on $I_1=[0,t]$}\quad \text{and}\quad \text{$r'\le-b\sqrt{r_{\max}^2-r^2}$ on $I_2=[t,y]$} \tag{7}
\end{equation}
-- or, equivalently, as
\begin{equation*}
\frac d{dx}\sin^{-1}\frac{r(x)}{r_{\max}}\ge b\text{ for }x\in(0,t),\quad
\frac d{dx}\sin^{-1}\frac{r(x)}{r_{\max}}\le -b\text{ for }x\in(t,y).
\end{equation*}
Integrating these inequalities over the corresponding intervals and recalling that $r(t)=r_{\max}$ and $r(y)=0$, we get
\begin{equation*}
\frac\pi2-\sin^{-1}\frac{r_0}{r_{\max}}\ge bt,\quad
-\frac\pi2\le-b(y-t);
\end{equation*}
eliminating $t$ from these two inequalities and recalling (5) (and the equality $r_0=f_0$), we finally get
\begin{equation*}
y\le\frac\pi b-\frac1b\,\sin^{-1}\frac{r_0}{r_{\max}}
\le\frac\pi b-\frac1b\,\sin^{-1}\frac1{\sqrt{1+\frac{g_{a,0}^2}{b^2f_0^2}}}=:y_1
\end{equation*}
-- in Case 1 (when $f_0=0$, let $y_1:=\frac\pi b$). Moreover, following the lines of the above reasoning, we see that if inequality (*) is replaced by the corresponding equality, then $y=y_1$ -- which shows that the upper bound $y_1$ on $y$ is exact in Case 1.
Case 2 is only simpler than Case 1. Here $r'<0$ on the entire interval $(0,y)$. Accordingly, here we need to consider only one branch of the inverse function, $p:=r^{-1}$, of the function $r$ on $[0,y]$. So, with $q:=p^2$, instead of (6) here we have
\begin{equation*}
g_{a,0}^2-{r'}^2=q(r_0)-q_j(r)\le-b^2(r_0^2-r^2)
\end{equation*}
for $r\in[0,r_0]$, that is, for $x\in[0,y]$.
So, in place of (7), here we have $r'\le-b\sqrt{r_*^2-r^2}$ on $[0,y]$, where $r_*$ is as in (5).
Thus, we similarly get
\begin{equation*}
y\le\frac1b\,\sin^{-1}\frac1{\sqrt{1+\frac{g_{a,0}^2}{b^2f_0^2}}}=:y_2
\end{equation*}
-- in Case 2. Moreover, following the lines of the above reasoning, we see that if inequality (*) is replaced by the corresponding equality, then $y=y_2$ -- which shows that the upper bound $y_2$ on $y$ is exact in Case 2.
One may also note that, if $g_{a,0}=0$ (that is, if $r'(0)=0$), then $y_1=y_2=\frac\pi{2b}$.
It is not hard to further generalize the result, by allowing $a=a(x)$ to depend on the argument $x$ of $f(x)$ and allowing $b=b(x,f(x))$ to depend on $x$ and $f(x)$ in such a way that $b_*:=\inf_{x,u} b(x,u)>0$.
Also, if in (*) we replace $\le0$ by $=H(x,f(x))$, where the function $H$ is such that $uH(x,u)\le0$ for all $u$, then it will follow that the resulting solution $f$ has infinitely many zeros and the distance between any consecutive strictly positive zeros will be bounded from above by $\pi/b_*$. Moreover, it will then follow that on each interval between consecutive strictly positive zeros the function $f$ will be concave if it is positive on such an interval, and $f$ will be convex on such an interval if it is negative on it.
