In addition to the above formulas, we can also use this very old formula, dating back to 1945 ( J. Ginsburg, Iterated exponentials, Scripta Math. 11 (1945), 340-353.):
$$f(x)=r+\sum_{n=1}^{\infty} \frac{\left(\ln a \right)^{n-1}\left(\ln \left(a^r \right)\right)^{nx}\left(1-r\right)^n B_n^x}{n!}$$B_n^{x-1}}{n!}$$
Where $B_n^x$ are the Bell numbers of x-th order and $r=\frac{W(-\log (a))}{\log (a)}$ ($W(x)$ is the Lambert function). Here: http://arxiv.org/abs/0812.4047 one can read about Bell numbers of higher orders.
The problem is that Bell numbers are only defined for integer order. We can easily generalize that to any real number by induction as follows:
$$A_0^x=1$$ $$A_{n+1}^x=\sum_{k=0}^{x-1} A_n^x\star A_n^k$$
And then $$B_n^x=A_{n-1}^{x+1}$$
where $f(n)\star g(n)$ is the binomial convolution as described by David Knuth:
$$f(n)\star g(n)=\sum_{k=0}^n \binom nkf(n-k)g(k)$$
To obtain the value for any real x, we can note that the right part in $A_{n+1}^x=\sum_{k=0}^{x-1} A_n^x\star A_n^k$ is a polynomial of x and k of degree n-1 and integer coefficients and we can take indefinite sum of it symbolically following the rule
$$\sum_x cx^n=\frac{B_{c+1}(x)}{c+1}$$
Where B_c(x) are the Bernoulli polynomials.
Unfortunately this method also works only for $a \le e^{1/e}$ in $f(x+1)=a^{f(x)}$.
Here is the plot of the function, for $a=\sqrt{2}$, obtained with this method and 5 terms:
