Edit: Re: whether or not there are R programs to do this - if it's normally distributed data with a pre-specified sample mean/covariance you're seeking to simulate, the mvrnorm function in R includes this option by setting empirical=TRUE. Here is a simple univariate example:
library(MASS)
x = mvrnorm(n = 10000, rep(0,1), 1, tol = 1e-6, empirical = TRUE)
mean(x)
[1] -5.793152e-18
var(x)
[,1]
[1,] 1
Original answer: To make your sample mean and variance exactly equal to a pre-specified value, you can appropriately shift and scale the variable. Specifically, if $X_1, X_2, ..., X_n$ is a sample, then the new variables
$$ Z_i = \sqrt{c_{1}} \left( \frac{X_i-\overline{X}}{s_{X}} \right) + c_{2} $$
where $\overline{X} = \frac{1}{n} \sum_{i=1}^{n} X_i$ is the sample mean and $ s^{2}_{X} = \frac{1}{n-1} \sum_{i=1}^{n} (X_i - \overline{X})^2$ is the sample variance are such that the sample mean of the $Z_{i}$'s is exactly $c_2$ and their sample variance is exactly $c_1$.
A similarly constructed example can restrict the range - shift/scale the data so that it lies within $(0,1)$ then multiply by $(b-a)$ and add $a$ to get it to lie in the interval $(a,b)$.