Derivation question of WKB method Quantum Mechanics (2nd Edition) by Bransden and Joachain contains the following passage:

Substituting (8.176) into (8.171), we obtain for $S(x)$ the equation
  $$-\frac{i\hbar}{2m}\frac{\mathrm{d}^2S(x)}{\mathrm{d}x^2} + \frac{1}{2}\biggl[\frac{\mathrm{d}S(x)}{\mathrm{d}x}\biggr]^2 + V(x) - E = 0.$$
  So far, no approximation has been made, this equation being strictly equivalent to the original Schrödinger equation (8.171). Unfortunately, equation 8.177) is a non-linear equation which is in fact more complicated than (8.171) itself! We must therefore try to solve (8.177) approximately. To this end, we first remark that if the potential is constant then $S(x) = \pm p_0 x$ (see (8.172)) and the first term on the left of (8.177) vanishes. Moreover, this term is proportional to $\hbar$, and hence vanishes in the classical limit ($\hbar\to 0$). This suggests that we treat $\hbar$ as a parameter of smallness and expand the function $S(x)$ in the power series
  $$S(x) = S_0(x) + \hbar S_1(x) + \frac{\hbar^2}{2}S_2(x) + \cdots\tag{8.178}$$
  Inserting the expansion (8.178) into (8.177) and equating to zero the coefficients of each power of $\hbar$ separately, we find the set of equations
  $$\begin{align}
\frac{1}{2m}\biggl[\frac{\mathrm{d}S_0(x)}{\mathrm{d}x}\biggr]^2 + V(x) - E &= 0\tag{8.179a} \\
\frac{\mathrm{d}S_0(x)}{\mathrm{d}x}\frac{\mathrm{d}S_1(x)}{\mathrm{d}x} - \frac{i}{2}\frac{\mathrm{d}^2 S_0(x)}{\mathrm{d}x^2} &= 0\tag{8.179b} \\
\frac{\mathrm{d}S_0(x)}{\mathrm{d}x}\frac{\mathrm{d}S_2(x)}{\mathrm{d}x} + \biggl[\frac{\mathrm{d}S_1(x)}{\mathrm{d}x}\biggr]^2 - i\frac{\mathrm{d}^2 S_1(x)}{\mathrm{d}x^2} &= 0\tag{8.179c}
\end{align}$$

My question is, why does each term have to be zero?
The equation (8.177) is equal to zero. After inserting (8.178) to (8.177).
It does not mean each term of order of $\hbar$ have to equal zero since $\hbar$ is constant, does it?
 A: In the semiclassical WKB approximation the Planck constant $\hbar$ is not a fixed number equal to its physical value $\approx 1.05 \times 10^{-34} Js$. Instead it is an indeterminate. The semiclassical WKB expansion is an asymptotic expansion in the limit $\hbar\to 0^+$.
A: I'm personally quite averse to treating $\hbar$ as "small parameter," since it has dimension and can be made arbitrarily big or small by suitable choice of units. Whenever this slight of hand is being done, what is really happening is that $\hbar$ is being treated as small compared to some other quantity that the author has elected not to identify. Thankfully, the WKB technique need not rely on any $\hbar\to0$ hocus pocus when done as an asymptotic expansion. Below I will work out the application of WKB to the Schrödinger equation. This post is based on Lecture 8 of Carl Bender's excellent guest lectures at Perimeter Institute (link to videos).
Starting with
$$
-\frac{\hbar^2}{2m}\psi''(x) = [E-V(x)]\psi(x) 
$$
let's multiply through by $-\frac{2m}{\hbar^2}$ and define $f(x) = \frac{2m}{\hbar^2}[V(x)-E]$. Then after letting $\psi=e^S$, our differential equation to tackle reads
$$
S''(x) + [S'(x)]^2 = f(x) \tag{$\star$}
$$
The physical motivation of this problem was to understand the wavefunction's behavior near the point $x_0$ where $E=V(x_0)$. Thus, $f\to0$ as $x\to x_0$. The WKB approximation is to say that near this turning point, $|S'(x)|^2\gg |S''(x)|$. For this post, let's take this as an inspired guess.
With the WKB approximation, we can replace the exact (but horrible) ODE ($\star$) with one that is asymptotically correct as $x\to x_0$
$$
[S'(x)]^2 \sim f(x)
$$
which has two solutions
$$
S(x) \sim \pm \int\sqrt{f(x)}dx$$
The sign "$\sim$" here stands for "is asymptotic to." Asymptotic relations follow many of the same rules as equality relations, but this post is not the place to get into that. Leave a comment and/or check out Prof. Bender's Lecture 7 for details if any manipulations below seem sketchy.
At this stage, we have the asymptotic behavior of $S(x)$ as $x\to x_0$, which means that we can write
$$
S(x) = \pm \int\sqrt{f(x)}dx + C(x) 
$$
and assert that
$$
|C(x)| \ll \left|\int\sqrt{f(x)}dx\right| \tag{$\star\star$}
$$
as $x\to x_0$. If we plug this new expression for $S(x)$ into ($\star$), we'll arrive at an ODE for $C(x)$:
$$
C'' + (C')^2 \pm 2\sqrt{f}C' \pm \frac{f'}{2\sqrt{f}} =0
$$
Disgusting! However, asymptotics will simplify life for us. From ($\star\star$) above, we can generate two more relationships
$$
|C'| \ll |\sqrt{f}| \tag{a}
$$
$$
|C''| \ll |f'/2\sqrt{f}| \tag{b}
$$
as $x\to x_0$. (a) lets us ignore the second term compared with the third, and (b) lets us ignore the first term compared with the fourth, leaving the relation (careful with the $\pm$ and $\mp$)
$$
C'(x) \sim -\frac{1}{4}\frac{f'(x)}{f(x)}
\implies
C(x) \sim -\frac{1}{4}\ln f(x)
\implies 
C(x) = -\frac{1}{4}\ln f(x) + D(x)
$$
with $|D(x)|\ll|C(x)|$ as $x\to x_0$.
Now, returning to our expression for $S(x)$, we can write
$$
S(x) ~ \pm\int\sqrt{f(x)}dx -\frac{1}{4}\ln f(x) + D(x)
$$
and exponentiating this, we get for the wavefunction
$$
\psi(x) = k\frac{e^{\pm \int\sqrt{f(x)}dx}}{\sqrt[4]{f(x)}}
$$
where $k = e^{D(x)}$ is a very slowly varying function that approaches a constant as $x\to x_0$.
A: 
It does not mean each term of order of $\hbar$ have to equal zero since $\hbar$ is constant, does it?

It's true that $\hbar$ is a known constant in the real world. But the theoretical model you're dealing with here is not quite the real world. It's a bit more general, and one of the ways in which it's more general is that it should work for any small value of $\hbar$, not just the actual value.
And when you have an equation of the form
$$C_0 + C_1 \hbar + C_2 \frac{\hbar^2}{2} + \cdots = 0$$
and you need it to be true for a whole range of values of $\hbar$, the only way for that to happen is that each of the coefficients $C_n$ is individually zero.
