Generating a nonrandom sequence which has a normal distributed density

Question

I need to create an algorithm in a computer program (Fortran90) which generates a sequence of $n$ (between $10$ and $10^6$) numbers $z$ that follow a normal distribution.

Restrictions:

1. Has to contain the number $1$.
2. Has to contain the number $z_i$, which is usually about $1500$, but has to be entered at runtime. That means that it has to be a variable in the algorithm and therefore it's not straightforward to parameterize.
3. It may not be done with random numbers.
4. The numbers do not have to be in order, even though it would be nice if they were.

To clarify: I need to find a function within a loop from $j=1$ to $n$ that calculates $z(j)$ in such a way that $z(k)=1$ ($k \epsilon \{1\dots n\}$), $z(l)=z_i$ ($l \epsilon \{1\dots n\}$), only a few numbers at the ''sides'' and more numbers ''in the middle''.
For example with $n=3$ and $z_i=1500$ the sequence would be $z=\{1,750,1500\}$.
For $n=5$, $z=\{1,a,750,b,1500\}$, where I have no idea how to calculate $a$ and $b$.

Background: I have to initialize a Lagrangian stochastic dispersion model (basically an air pollution model) with $n$ particles in height $z(n)$. Usually they are all started at a source in a given height, but for testing purposes I need them to start with a normal distributed density: Few at the bottom, few at the top and many in the middle. It would be easy to do this randomly distributed, but I need to ensure that the start is normal distributed even for a small amount of particles.

Disclaimer: I'm not sure whether this is the correct stackexchange, nomenclature in the title or tags.

EDIT: This is a general version of my solution based on majlon's answer here:

zi=upper value (the 1500 above)
zl=lower value (the 1 above)
n=number of elements in z
z=desired values as described above
probit=inverse cumulative distribution function 
       associated with the standard noral distribution
--------
mean=(zi-zl)/2
std=(zl-mean)/probit(1/n/2.)
do i=1,n
   z(i)=mean+probit( (2.*i-1)/2./n) )*std
end do

Answer 1

You need to utilize the percentiles of the inverse CDF of the gaussian-distribution described as Wikipedia. Divide the quantiles over the amount of numbers wanted and then adapt $\sigma$ so that the outer segments align with 1 and Z ($\mu$ will be (Z-1)/2).

So if n=3, that means that each segment is 33%. That means that 1 and Z should be at in the middle of (0-33) and (67-100) segments respectively. This gives $\mu + \sigma\Phi^{-1}(0.33/2)=1$. If you solve $\sigma$ you can find a gaussian distribution that would fit and then use the other corresponding percentiles to get your other values. I am unsure how you can get the error function in Fortran90 but it is doable as long as you have an approximation of e.

I hope that gets you on your way.