The function you want grows too fast to be interpolated by usual method, but there exists an iterative solution with Cauchy integrals by Dmitry Kouznetsov
If you relax the condition so to find a solution for $f(x+1)=a^{f(x)}$ such that $$a \le e^{1/e} $$ then there are multiple expressions for your function:
$$f(x)=\sum_{m=0}^{\infty} \binom xm \sum_{k=0}^m \binom mk (-1)^{m-k}\exp_a^{[k]}(1)$$
$$f(x)=\lim_{n\to\infty}\binom xn\sum_{k=0}^n\frac{x-n}{x-k}\binom nk(-1)^{n-k}\exp_a^{[k]}(1)$$
$$f(x)=\lim_{n\to\infty}\frac{\sum_{k=0}^{2n} \frac{(-1)^k \exp_a^{[k]}(1)}{(x-k)k!(2n-k)!}}{\sum_{k=0}^{2n} \frac{(-1)^k }{(x-k) k!(2n-k)!}}$$
$$f(x)=\lim_{n\to\infty} \log_a^{[n]}\left(\left(1-\left(\ln \left(\frac{W(-\ln a)}{-\ln a}\right)\right)^x\right)\frac{W(-\ln a)}{-\ln a}+\ln \left(\frac{W(-\ln a)}{-\ln a}\right)\exp_a^{[n]}(1)\right)$$
Always here the number in square brackets designates n-th iteration and $W(x)$ is the Lambert's function.
There is also expression for inverse function:
$$ f^{[-1]}(x)=\left(\lim_{n\to\inftyf^{[-1]}(x)=\lim_{n\to\infty} \frac{\ln \left(\frac{\frac{W(-\ln a )}{\ln a}+\exp_a^{[n]}(x)}{\frac{W(-\ln a)}{\ln a}+\exp_a^{[n]}(1)}\right)}{\ln \ln \left(\frac{W(-\ln a)}{-\ln a}\right)}\right)$$a}\right)}$$
