A very similar question was asked even earlier in Atkin–Lehner operator for GL(3)?. I was looking for the same generalisations of the Atkin-Lehner operators for $GL(n)$ and I managed to find something useful recently. Though the question is old, I think it could help other people like me to give an answer here.
In short, the operator you are looking for is given by $$ Wf(z)=f(\text{diag}(1,1,N) \cdot z^{-T}).
$$
This is the operator that gives you the functional equation and it generalises all the other properties you would want from a Fricke involution.
I was inspired in this definition by the dual form $\tilde{f}(z)=f(wz^{-T}w^{-1})$, where $w$ is the long Weyl element. For $SL_n(\mathbb{Z})$, this gives the functional equation, but it is no longer an isomorphism when passing to the congruence subgroups (it does not leave the group invariant). It turns out that we need to use the dual form because the normaliser of $\Gamma_0(N)$ is trivial in higher rank and cannot give rise to any interesting operators.
In fact, you can try to generalise this approach to get other Atkin-Lehner operators, like the ones corresponding to certain divisors of $N$. Unfortunately, the fact that there is a shortage of Atkin-Lehner operators for powerful levels is a phenomenon that blows up when going to higher rank and, at least according to my definitions, there are no other operators apart from the Fricke involution for $n>2$.
I wrote a note giving more details and proofs for my claims. You can access it on my website at https://www.math.uni-bonn.de/people/toma/notes/AL-operators.pdf. It's still rough around the edges, but I hope it gives a satisfying answer, and I'd be interested in comments or suggestions.