A procedural solution for $n$ tuple convolution is given here.
Q1: The two gamma distribution convolution (GDC) in closed form is only given here
$\mathrm{G}\mathrm{D}\mathrm{C}\left(\mathrm{a}\kern0.1em ,\mathrm{b}\kern0.1em ,\alpha, \beta; \tau \right)=\left\{\begin{array}{cc}\hfill \frac{{\mathrm{b}}^{\mathrm{a}}{\beta}^{\alpha }}{\Gamma \left(\mathrm{a}+\alpha \right)}{e}^{-\mathrm{b}\tau }{\tau^{\mathrm{a}+\alpha-1}}{}_1F_1\left[\alpha, \mathrm{a}+\alpha, \left(\mathrm{b}-\beta \right)\tau \right],\hfill & \hfill \tau >0\hfill \\ {}\hfill \kern2em 0\kern6.6em ,\hfill \kern5.4em \tau \kern0.30em \le \kern0.30em 0\hfill \end{array}\right.,$
which is a density function consisting of a gamma variate multiplied by $_1 F _1(A; B; Z)$, where the latter is a confluent hypergeometric function of the first kind, and see end note below. For $\text{b} = β$, this equation reduces to the well known result
$$\text{iff b}=\beta\text{;}\;{\mathrm{GD}}\left(\mathrm{a},\mathrm{b};\;\tau \right)\otimes {\mathrm{GD}}\left(\alpha, \text{b}; \tau \right)=\text{GD}(\text{a}+α,\text{b} ;τ)$$
An even more general solution with weights is given in Eq. (2) here. Furthermore, that reference lists a non-closed form solution for the sum of $n$ weighted gamma distributions.
Q2: No applications for the GDC have been posted on CV prior to this. Most recently, the GDC has been used in medicine for the first time to model radioactive tracer activity in time in the thyroid gland. GDC models have also been applied for ecological water storage I/O, waiting times in queuing theory, and in the evaluation of aggregate economic risk of portfolios.
End note: The fast computation of the confluent hypergeometric function of the first kind uses the Euler's integral identity ${}_1F_1\left(A;B;Z\right)=\frac{\Gamma \left(B\right)}{\Gamma \left(B-A\right)\Gamma \left(\mathrm{A}\right)}{\displaystyle {\int}_0^1{\mathrm{e}}^{Z\;u}{u}^{A-1}{\left(1-u\right)}^{B-A-1}du}$.
Best Answer
Let $X$ and $Y$ be independent random variables with densities $f, g$ respectively. Then the convolution of $f$ and $g$ $$ (f * g )(t) \, \stackrel{\mathrm{def}}{=}\ \int_{-\infty}^\infty f(\tau)\, g(t - \tau) \, d\tau $$ is the density function of $X+Y$. For this to make dimensional meaning, we must assume that $X$ and $Y$ are measured with the same units of measurement, it could for instance be m/s (meters per second). Then the unit of measurement of the probability density function is probability per (m/s). Since probability is a pure number we can write this as 1/(m/s). Let us make this more general by writing u for whatever common unit of measurement of $X$ and $Y$, then the densities $f,g$ has units 1/u.
We can represent the convolution integral above with a Riemann sum: $$ = \sum_i f(\tau_i) g(t-\tau_i) (\tau_i -\tau_{i-1}) $$ The variable $\tau$ above clearly has unit u. So the unit of measurement of each term in the Riemann sum has unit $$ \frac{1}{\text{u}}\cdot \frac{1}{\text{u}}\cdot\text{u} = \frac{1}{\text{u}} $$ showing that the density function obtained by convolving $f$ and $g$ has the same unit as $f$ and $g$.