Neyman Allocation (or modifications of NA) is often used in practice. Yes, you are right, we never know $S_h$ when doing the calculation of sample allocation. But we can estimate $S_h$ or use some approximation of $S_h$.
Assume $S_h(y)$ is computed for a variable $\textbf{y}$. $S_h$ can be estimated from the previous survey, if $\textbf{y}$ was observed in the previous survey.
There could be another variable $\textbf{z}$ which is correlated with $\textbf{y}$. $\textbf{z}$ could be available from other survey or some auxiliary data source (register, census). Then you can use $S_h(z)$ or $s_h(z)$ (an estimate) as approximate for $S_h(y)$.
It could be possible to guess $S_h(y)$, for example in case of binary $\textbf{y}$.
Keep in mind - you will never achieve optimal allocation in practice, so allocation close to the optimal could be good enough.
The quoted formula is not quite right. Let's derive the correct one.
Since the population mean (or any other constant) may be subtracted from every value in a population $S$ without changing the variance of the population or of any sample thereof, we might as well assume the population mean is zero. Letting the values in the population be $\{x_i\, \vert\, i\in S\}$, this implies
$$0 = \sum_{i\in S} x_i.$$
Squaring both sides maintains the equality, giving
$$0 = \sum_{i,j\in S}x_ix_j = \sum_{i\in S}x_i^2 + \sum_{i \ne j \in S} x_ix_j,$$
whence
$$\sum_{i\ne j \in S} x_ix_j = -\sum_{i\in S} x_i^2.$$
This key result will be employed later.
Let $S$ have $N$ elements. Because its mean is zero, its variance is the average squared value:
$$s^2 = \frac{1}{N}\sum_{i\in S}x_i^2.$$
(Please note that there can be no dispute about the denominator of $N$; in particular, it definitely is not $N-1$: this is a population variance, not an estimator.)
To find the variance of the sample distribution of the mean, consider all possible $n$-element samples. Each corresponds to an $n$-subset $A\subset S$ and has mean
$$\frac{1}{n}\sum_{i\in A} x_i.$$
Since the mean of all the sample means equals the mean of $S$, which is zero, the variance of these $\binom{N}{n}$ sample means is the average of their squares:
$$s_n^2 = \frac{1}{\binom{N}{n}} \sum_{A\subset S}\left(\frac{1}{n}\sum_{i\in A}x_i\right)^2 = \frac{1}{n^2\binom{N}{n}} \sum_{A\subset S}\sum_{i,j\in A}x_ix_j \\= \frac{1}{n^2\binom{N}{n}} \sum_{A\subset S}\left(\sum_{i\in A}x_i^2 + \sum_{i\ne j\in A}x_ix_j\right) .$$
(Once again, $\binom{N}{n}$, not $\binom{N}{n}-1$, is the correct denominator: this is the variance of a collection of $\binom{N}{n}$ numbers, not an estimator of anything.)
Fix, for a moment, any particular index $i$. The value $x_i$ will appear in $\binom{N-1}{n-1}$ samples, because each such sample supplements $x_i$ with $n-1$ more elements of $S$ out of the $N-1$ remaining elements (sampling is without replacement, remember). Its contribution to the right hand side therefore equals $\binom{N-1}{n-1}x_i^2$.
Also fixing an index $j\ne i$, similar reasoning shows the product $x_ix_j$ appears in $\binom{N-2}{n-2}$ samples, thereby contributing $\binom{N-1}{n-1}x_ix_j$ to the right hand side. Therefore, upon summing over all such $i$ and $j$ in $S$,
$$s_n^2 = \frac{1}{n^2\binom{N}{n}} \left(\binom{N-1}{n-1}\sum_{i\in S}x_i^2 + \binom{N-2}{n-2}\sum_{i\ne j\in S}x_ix_j\right).$$
Plug the first result into that last sum:
$$s_n^2 = \frac{1}{n^2\binom{N}{n}} \left(\binom{N-1}{n-1}\sum_{i\in S}x_i^2 + \binom{N-2}{n-2}\left(-\sum_{i\in S}x_i^2\right)\right).$$
It is now straightforward to relate this to the variance of $S$, because $\sum_{i\in S}x_i^2 = Ns^2$:
$$s_n^2 = \frac{1}{n^2\binom{N}{n}} \left(\binom{N-1}{n-1} - \binom{N-2}{n-2}\right)\left(Ns^2\right) = \frac{s^2}{n}\left(1 - \frac{n-1}{N-1}\right).$$
Thus the sampling variance for sampling with replacment, $\frac{s^2}{n}$, is multiplied by $1 - \frac{n-1}{N-1}$ to obtain the sampling variance for sampling without replacement, $s_n^2$. Accordingly, the multiplicative adjustment for the sampling standard deviation is its square root, $\sqrt{1- \frac{n-1}{N-1}}$. This differs from the quoted formula, which uses $\sqrt{1 - \frac{n}{N}}$.
Two simple checks can give us some comfort concerning the correctness of this result. First, the sample variance of means of samples of size $n=1$, $s_1^2$, obviously equals the population variance $s^2$. The correct formula states
$$s_1^2 = \frac{s^2}{1}\left(1 - \frac{1-1}{N-1}\right) = s^2,$$
as it should. Unfortunately, the quoted formula asserts that $s_1^2 = s^2(\frac{1}{1} - \frac{1}{N})$ which obviously cannot be right. Second, the sample variance of the means of samples of size $n=N$ is zero, because there is no variation, and indeed both formulas give $0$ in this case.
Best Answer
I'm assuming you are familiar with the basic theory/mathematics underlying the formulas to obtain the variance when working with simple random samples with a finite population. When working with a finite population, then under simple random sampling (without replacement), the variance of the mean is given by:
$V(\bar{Y}_{SRS})=\frac{1-f}{n}S^{2}$
Under stratified sampling the mean is given by: $\sum_{h=1}^{L}\frac{N_{h}}{N}\bar{Y}_{h}$, which is just a weighted sum of the individual stratum means. So we can compute the variance of the stratified mean as follows, since, under stratified random sampling, sampling is performed independently in each stratum and therefore the variance of the sum of the stratum means is the sum of the variances of the stratum means (i.e. the covariance term vanishes):
\begin{eqnarray*} V\left[\sum_{h=1}^{L}\frac{N_{h}}{N}\bar{Y}_{h}\right] & = & \sum_{h=1}^{L}V\left(\frac{N_{h}}{N}\bar{Y}_{h}\right)\\ & = & \sum_{h=1}^{L}\left(\frac{N_{h}}{N}\right)^{2}V\left(\bar{Y}_{h}\right)\\ & = & \sum_{h=1}^{L}\left(\frac{N_{h}}{N}\right)^{2}\frac{1-f_{h}}{n_{h}}S_{h}^{2}\\ & = & \sum_{h=1}^{L}\left(\frac{N_{h}}{N}\right)^{2}\left(1-f_{h}\right)\frac{S_{h}^{2}}{n_{h}}\\ & = & \sum_{h=1}^{L}\left(\frac{N_{h}}{N}\right)^{2}\left(1-\frac{n_{h}}{N_{h}}\right)\frac{S_{h}^{2}}{n_{h}}\\ & = & \sum_{h=1}^{L}\left(1-\frac{n_{h}}{N_{h}}\right)\left(\frac{N_{h}}{N}\right)^{2}\frac{S_{h}^{2}}{n_{h}}\blacksquare \end{eqnarray*}