I suspected something deeper was going on. Indeed, if we just substitute $L$ into the recursions, we get trivial equations. So clearly we can't work our way out that easily!
I've often seen situations where invariants come in useful. Given a sequence $x_n$, a function $f$ such that $f(x_n) = f(x_{n+1})$ for all $n$ is called an invariant of the sequence.
The key thing, is that if $f$ is continuous then $f(x_0) = f(x_1)=\ldots = f(\lim_{n \to \infty} x_n)$, and therefore we may be able to use the nature of preimages of $f$ to calculate our answer.
The first thing I did was write Python code to simulate the situation. I just took some arbitrary $a,b,c$ here.
from math import *
def checkone() :
a = 3
b = 7
c = sqrt(61)
for i in range(100) :
d = a
e = b
f = c
a = sqrt(d*(f+e-d))
b = sqrt(e*(f+d-e))
c = sqrt(f*(d+e-f))
print(str([a,b,c]) + ","+str(a**4+b**4+c**4))
checkone()
where I used the $a^4+b^4+c^4$ term to predict that something was going on with the fourth powers of the sequence. I noticed this fourth-power pattern when I saw that we would need to play with quadratics of $x_n^2$, which due to the nature of the RHS (the fact that $x_{n+1}^2 = sth. x^2 + sth. x_ny_n $ etc.) meant that I'd have to go a square further.
Call $(L(a,b,c),L(a,b,c),L(a,b,c))$ as the limit obtained from the above. Following some playing around, I first came to the conclusion that :
$$
3L^4(a,b,\sqrt{a^2+b^2}) = 4a^2b^2
$$
(note : I actually used the homogeniety of the relation and realized that it was enough to take $a=1$, so this simplified things) and then by perturbation to the conclusion:
$$
3L^4(a,b,\sqrt{a^2+b^2-c}) = 4a^2b^2-c^2
$$
This then tells me what the limit should be , following a change of subject:
$$
3L^4(a,b,x) = 2(a^2b^2+b^2x^2+x^2a^2) - a^4-b^4-x^4 \\
$$
Now comes the wonderful idea : That quantity, right there on the RHS, is actually invariant.
Suppose that $x_n,y_n,z_n$ are defined as given. Then :
$$
2(x_{n+1}^2y_{n+1}^2+x_{n+1}^2z_{n+1}^2+y_{n+1}^2z_{n+1}^2) - x_{n+1}^4-y_{n+1}^4-z_{n+1}^4\\ = 2(x_n^2y_n^2+y_n^2z_n^2+z_n^2x_n^2)-x_n^4-y_n^4-z_n^4
$$
This is a non-routine but fairly ok check to make. So you substitute the recursion where required and make sure things simplify.
Once we make this check, we note that we can, obviously by induction, obtain the statement :
$$
2(x_{n}^2y_{n}^2+x_{n}^2z_{n}^2+y_{n}^2z_{n}^2) - x_{n}^4-y_{n}^4-z_{n}^4\\ = 2(x_0^2y_0^2+y_0^2z_0^2+z_0^2x_0^2)-x_0^4-y_0^4-z_0^4
$$
for all $n \geq 1$. Now, the splendid work of the OP has proved that the sequence converges to a limit $(\alpha,\alpha,\alpha)$.
But now, we may take the limit on both sides of the recursion. Indeed, note that if we do so then we get :
$$
3L(x_0,y_0,z_0)^4 = 2(x_0^2y_0^2+y_0^2z_0^2+z_0^2x_0^2)-x_0^4-y_0^4-z_0^4
$$
and therefore :
$$
\boxed{
L(x_0,y_0,z_0) = \sqrt[4]{\frac{2(x_0^2y_0^2+y_0^2z_0^2+z_0^2x_0^2)-x_0^4-y_0^4-z_0^4}{3}}
}
$$
which gives the correct answer for the case above.
So find invariants, is the lesson. Don't ask me how I'd do this without computing power, though!
Best Answer
The following two descriptions of the upper and lower limit might turn out to be more useful here:
$x \in \mathbb R \cup \{ \pm \infty \}$ is a limit point of $(x_n)$ if and only if there exists a subsequence $(x_{n_k})$ with $x_{n_k} \to x$;
$\liminf x_n$ is the infimum of all the limit points of $(x_n)$; in other words, there exists a subsequence $(x_{n_k})$ with $x_{n_k} \to x$ and if there exists $y$ and a subsequence $(x_{m_l})$ with $x_{m_l} \to y$ then $x \le y$;
a similar characterization for $\limsup$;
if you work with nets instead of sequences, then you'll have to use subnets instead of subsequences; everything else remains unchanged.
1) If $\liminf z_n = -1$, there exists a subsequence $(z_{n_k})$ with $z_{n_k} \to -1$. It follows immediately that the subsequence $(x_{n_k} + y_{n_k} z_{n_k})$ of the sequence $(x_n + y_n z_n)$ converges to $-\infty + (-1) \cdot \infty = -\infty$, showing that $-\infty$ is a limit point of $(x_n + y_n z_n)$. Since, by its very nature, $-\infty$ is the smallest possible element of $\mathbb R \cup \{ \pm \infty \}$, every other limit point of $(x_n + y_n z_n)$ must be greater than or equal to it, therefore $-\infty$ is the infimum of all the limit points of $(x_n + y_n z_n)$, therefore is its $\liminf$.
2) Exactly as above, with $-\infty + (-1) \cdot \infty$ getting replaced by $-\infty + (-1) \cdot 0$, which is still $-\infty$.
Notice that, so far, we haven't needed the monotonicity of $(x_n)$ and $(y_n)$, but only their limits. Notice also that we haven't used that $-1$ is the $\liminf$ of $(z_n)$, but only that it is a limit point.
3) Again, with the same construction as above, the subsequence $(x_{n_k} + y_{n_k} z_{n_k})$ tends to $0 + (-1) \cdot 0 = 0$, so $0$ is a limit point of $(x_n + y_n z_n)$. If $r < 0$ is any other limit point, then there exists a subsequence $q_{m_l}$ of $(x_n + y_n z_n)$ with $q_{m_l} \to r$. Consider then the subsequences $(x_{m_l})$, $(y_{m_l})$ and $(z_{m_l})$. Since $(x_n)$ is convergent (to $0$), so will be $(x_{m_l})$, therefore $y_{m_l} z_{m_l} \to r$.
If $(y_n)$ is constant $0$ from some $N$ onwards, then $x_n + y_n z_n = z_n$ from that $N$ onwards, and the result is trivial. If $(y_n)$ is not eventually $0$, then from $(y_{m_l})$ I may extract yet another subsequence $(y_{m_{l_i}})$ such that $y_{m_{l_i}} \ne 0$ for all $i$. Since $(y_n)$ decreases to $0$, it follows that $y_n \ge 0$, therefore $y_{m_{l_i}} > 0$.
Since $y_{m_{l_i}} z_{m_{l_i}} \to r$, we deduce that for every $\varepsilon > 0$ there exists $I \in \mathbb N$ such that for $i \ge I$ we have $| y_{m_{l_i}} z_{m_{l_i}} - r | < \varepsilon$. Taking $\varepsilon = - \frac r 2$, this implies that for $i \ge I$ we have $y_{m_{l_i}} z_{m_{l_i}} - r < \frac r 2$, which implies that $z_{m_{l_i}} < \frac r {2 y_{m_{l_i}}} \to -\infty$ (because $r<0$), which implies that $zy_{m_{l_i}} \to -\infty$, which means that $-\infty$ is a limit point of $(z_n)$. But this is a contradiction, because $-1$ is its $\liminf$, which is the infimum of the limit points. Therefore, $r \ge 0$, i.e. every other limit point of $(x_n + y_n z_n)$ is $\ge 0$, which (together with the already proven fact that $0$ is a limit point) makes $0$ its $\liminf$.
Notice that, unlike for (1) and (2), this time we have used the monotonicity of $(y_n)$ in an essential way, but (again) not the one of $(x_n)$. Notice also that we have used the fact that $-1$ is the $\liminf$ of $(z_n)$ in an essential way (unlike for (1) and (2)).