I'm trying to solve the following equation with FVM
$$u_t + f(u)_x = g(u)$$
where $g$ is some smooth function of $u$ and $f(u) = \frac{u^2}{2}$. This is really similar to Burger's equation, except for the source term. My main problem is that I don't know how to derive correctly the numerical scheme since I have a source term here.
More precisely, taking a space-time cell:
$\int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} u(t_{n+1},x) - u(t_{n},x) + \int_{t_n}^{t_{n+ 1}} f(u(x_{j+\frac{1}{2}}),t) - f(u(x_{j-\frac{1}{2}},t))dt = \int_{t_n}^{t_{n+1}} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} g(u)du$
Now I define as usual \begin{align} \bar{u}_j^n = \frac{1}{\Delta x} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} u(x,t_n)dx \end{align} and
\begin{align} F_{j+ \frac{1}{2}^n} = \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} f(u(x_{j+\frac{1}{2}},t))dt \end{align}
and hence the numerical scheme writes
$$\bar{u}_{j+1}^n = \bar{u}_j^n - \frac{\Delta t}{\Delta x} (F_{j+ \frac{1}{2}^n} - F_{j-\frac{1}{2}}^n) + g(\bar{u}_j^n)$$
where I replaced the average of $g(u)$ with its evaluation on the average, which is a second order approximation.
My question: in order to determine the numerical fluxes I usually solve Riemann problems at the interfaces $x=x_{j+\frac{1}{2}}$, like
$u_t + f(u)_x = 0$ with initial data $u(x,t_n)$ given by \begin{align} \bar{u}_{j-1}^n , \quad x< x_{j-\frac{1}{2}}\\ \bar{u}_j^n, \quad x> x< x_{j+\frac{1}{2}} \end{align}
My problem is that now I should also include the source term in the Riemann problem, so how can I handle it? I need to understand how to find an expression for the numerical fluxes in this case