Analogous to the Lambert W function, define $H(x)$ be the inverse function of $x e^{x e^x}$, i.e. $$
H(x) e^{H(x)e^{H(x)}} = x
$$
for $x\in\mathbb{R}$. This is well defined for all real $x$, since $xe^{xe^x}$ is monotonic. The function you are looking for can be found by observing $$
e^{H(\log s)e^{H(\log s)*e^{H(\log s)}}} = s
$$
so $e^{H(\log s)}$ is the 3rd super-root of $s$, analogous to how $e^{W(\log s)}$ is the super square root of $s$.
Just as $W(t)$ can be written as the integral of an elementary function$$
W(t) = \frac{t}{\pi} \int_0^\pi\frac{(1-x\cot x)^2 + x^2}{t + x\cot x e^{-x\cot x}}dx
$$
we will find an analogous formula for $H$:
$$
H(t) = \frac{t}{\pi}\int_{-\infty}^\infty \Im\left(\frac{1 + ((x+i\frac\pi4)^2 + x+i\frac\pi4)e^{x+i\frac\pi4}}{t - (x+i\frac\pi4)e^{(x+i\frac\pi4)e^{x+i\frac\pi4}}}\right) dx
$$
(See WolframAlpha if you want to know what the integrand looks like in closed form without $i$ or imaginary parts.) See below for a proof. This is not particularly nice, but it means that $x^{x^x}=y$ can be solved for $x$ using only elementary functions and integration. As a numerical check, suppose we want to solve $x^{x^x} = 2$. Then we would take $x = e^{H(\log 2)}$. Using WolframAlpha, we compute that $H(\log 2) \approx 0.389799$. Exponentiating gives $x=1.47668$. This is not far off:$$
1.47668^{1.47668^{1.47668}} = 1.99998...
$$
Just to show that it also works for at least some values between 0 and 1, lets also solve $x^{x^x} = \frac12$. We find $H(-\log 2) = -1.00137$, hence we should have $x \approx 0.36736$ is the solution. Indeed:$$
0.36736^{0.36736^{0.36736}} = 0.499984
$$
We prove the integral formula for $H$ using contour integration: Fix $t$. Observe that the function $$
f(x) = \frac{1 + (x+x^2)e^x}{xe^{xe^x} - t}
$$
has exactly one pole on the real line, at the point $x=H(t)$. By noting that $\frac{d}{dx} xe^{xe^x} = (1+(x+x^2)e^x)e^{xe^x}$, we find that the residue at this pole is $\frac1{e^{{H(t)}e^{H(t)}}} = \frac{H(t)}t$. Hence for any counterclockwise oriented curve in $\mathbb{C}$ that surrounds $H(t)$ and includes no other solutions to $ze^{ze^z}=t$ in its interior, we have $$
\frac{1}{2\pi i}\oint_\gamma f(x)dx = \frac{H(t)}{t}
$$
Take $\gamma$ to be a rectangular contour:
Note that the integral along the top minus the integral along the bottom is equal to $$
2\int_{-A}^A \Im(-f(x+i\frac\pi4))dx
$$
(negative because we are integrating right to left on the top of the box). Thus we need to show that the integral over the left and right sides of the box goes to 0, and that $f$ has no more poles inside the region. To see that $f$ goes to 0 on the left and right sides, observe that as $x$ goes to infinity, if $|\Im x|<\pi/2$, then the denominator goes to infinity doubly exponentially, making the whole fraction go to 0 rapidly. When $x$ goes to negative infinity, the denominator goes to infinity linearly, while the numerator is bounded. Since both the right and left legs of the contour have constant length, if $f$ goes to 0 on them as $A$ goes to infinity, these do not contribute to the integral in the limit.
Thus it remains to show that $f$ has no other poles in the interior of the contour. Letting $z=x+i y$, we compute $$
\arg {z e^{z e^z}} = \arg z + \Im(z e^z) = \arg z + e^x x \sin(y) + e^x y \cos(y)
$$
For $z$ in the subregion defined by $\Re z>-1$, $\Im z \in(0,\pi/4)$, all terms are positive, so $z e^{z e^z}$ is guaranteed to be nonreal. Note also that for $\Re z < -1$, we have $$
|z e^{z e^z}| > e^{-|z e^z|} = e^{-|z| e^x} \ge e^{-|x|e^x -\frac\pi4 e^x} \ge c
$$
we $c=\min\limits_{x<-1} e^{-|x|e^x -\frac\pi4 e^x} = \exp\left(\frac{-4-\pi}{4e}\right)\approx 0.5185$. Hence if $|t|<0.5185$, $f$ cannot possibly have any poles in the region of integration, so we can solve $$
x^{x^x} = y
$$
using this method at least for $|\log x|< c$. It should also work in a wider range, by analytic continuation in the largest interval where the integral converges it ought to give correct values for $H$, since the integrand is an analytic function of $t$. Numerically it seems to fail for $x=\frac18$.
Let's consider solutions in closed form.
$$e^x\ln(x)=a$$
1.) No algebraic solutions except $1$ for algebraic $a$
According to Schanuel's conjecture, we have $\forall\ 0,1\neq x\in\overline{\mathbb{Q}}\colon$ $\text{trdeg}_\mathbb{Q}\mathbb{Q}\left(x,\ln(x),e^x\right)=2$. That means: $\ln(a_0),e^{a_0}$ are $\overline{\mathbb{Q}}$-algebraically independent for each algebraic $a_0\neq 0,1$. That means, your equation cannot have algebraic solutions except $1$ for algebraic $a$.
2.) No general elementary solutions
According to the theorem of [Ritt 1925], that is proved also in [Risch 1979], the function $f\colon x\mapsto e^x\ln(x)$ seems not to have partial inverses that are elementary functions for non-discrete domains of $f$.
But that doesn't say if solutions in the elementary numbers are possible for single $a$.
[Risch 1979] Risch, R. H.: Algebraic Properties of the Elementary Functions of Analysis. Amer. J. Math. 101 (1979) (4) 743-759
[Ritt 1925] Ritt, J. F.: Elementary functions and their inverses. Trans. Amer. Math. Soc. 27 (1925) (1) 68-90
3.) Solutions in terms of Hyper Lambert W
$$e^x\ln(x)=a$$
$$x\to e^z: e^{e^z}(z+2k\pi i)=a\ \ \ \ \ (k\in\mathbb{Z})$$
$$-\pi<\text{Im}(z)\le\pi:$$
$$e^{e^z}z=a$$
This equation can be solved by Hyper Lambert W:
Galidakis, I. N.: On solving the p-th complex auxiliary equation $f^{(p)}(z)=z$. Complex Variables 50 (2005) (13) 977-997
Galidakis, I. N.: On some applications of the generalized hyper-Lambert functions. Complex Variables and Elliptic Equations 52 (2007) (12) 1101-1119
Best Answer
As @Zoe Allen commented, squaring and letting $z=\cos(x)$, we have to solve for $z$ the equation $$y^2=(1-z^2)\,e^{2z}\qquad \implies\qquad e^{-2z}=\frac {1-z^2} {y^2}$$ which has a solution in terms of the generalized Lambert function (have a look at equation $(4)$).