You can do this in R with optim, e.g. minimizing the sum of squared errors for the mean and the 95th percentile:

obs = c(28.744,385.714,20.595,99.350,31.864,77.713,       
264.408,21.204,31.937,0.900,18.762,173.276,23.707)
tail = 500.912

error = function(x)(qgamma(0.95, shape=x[1], scale=x[2]) - tail)^2 + (x[1]*x[2]-mean(obs))^2

optim(c(1,1), error)

This uses the formula for the mean and gets the desired mean and percentile exactly with a shape parameter of 0.147 and a scale parameter of 616.6.

You can fit a gamma distribution in R with optim, e.g. by minimizing the sum of squared errors for the mean and the 95th percentile, which actually sets both errors to zero:

obs = c(28.744,385.714,20.595,99.350,31.864,77.713,       
264.408,21.204,31.937,0.900,18.762,173.276,23.707)
tail = 500.912

error = function(x)(qgamma(0.95, shape=x[1], scale=x[2]) - tail)^2 + (x[1]*x[2]-mean(obs))^2

params = optim(c(1,1), error)$par
params

The code uses the formula for the gamma’s mean and gets the desired mean and percentile exactly with a shape parameter of 0.147 and a scale parameter of 616.6.

You can then check the quantile with

qgamma(0.95, shape=params[1], scale=params[2])

You can also check the log-likelihood of the result with

sum(log(dgamma(shape = params[1], scale = params[2], obs)))

and compare it with the log-likelihood from other two-parameter distributions.