update: an update on the Pari/GP routine, few textual improvements
I've got a sudden enlightenment which seems to solve the problem completely.
The key here is, that on one hand the fixed- and periodic points of $\exp(z)$ are all repelling (this was also mentioned/proved in an article of Shen/Rempe-Gillen, pg 1, see reference below), and on the other hand the one single fixed point of $\log(z)$ is nicely attracting, but it is only a singular occurence.
$1$-periodic points/fixed points
Well, further fixed points can be found by using branches of the $\log()$, say
$$ \text{lb}(z,k) = \log(z) + k \cdot C \qquad \text{where } C = i 2 \pi$$
We may then have for some $p_1 = p_1' + C $ the equality
$$ \exp(p_1) = p_1 = p_1' + C = \exp (p_1'+ C) $$
and on the other hand the inverse relation
$$ \text{lb} (p_1,1) = p_1'+ C = \text{lb}(p_1'+C,1)$$
Noticing this we can exploit the attractive property of the (iterated) logarithmizing to converge from some initial value to the desired fixpoint, say from $z=1+I$ to some $p_1$ .
Example:
z = 1+I | %681 = 1 + I
z=lb(z,1) | %682 = 0.34657359 + 7.0685835*I
z=lb(z,1) | %683 = 1.9568606 + 7.8049907*I
z=lb(z,1) | %684 = 2.0852451 + 7.6083262*I
z=lb(z,1) | %685 = 2.0654578 + 7.5864764*I
From where we can now use the Newton-iteration to get even quadratic convergence rate and then find the arbitrary well approximated value of the desired fixed point $\exp(z)=z$
z=newton(z) | %686 = 2.0622777 + 7.5886312*I
exp(z)-z | %687 = 6.6538205 E-211 - 1.9010916 E-211*I
Using the $k$'th multiple of the constant $C=i2\pi$ we can now find (and index) the (infinitely many) complex fixed points ($1$-periodic):
$$ p_1 = \text{lb}(p_1,k) \to \exp(p_1)=p_1$$
Another example:
p_1 = 1+I | %688 = 1 + I
p_1=lb(p_1,3) | %689 = 0.34657359 + 19.634954*I
p_1=lb(p_1,3) | %690 = 2.9774671 + 20.402703*I
p_1=lb(p_1,3) | %691 = 3.0262041 + 20.275440*I
p_1= newton(p_1) | %692 = 3.0202397 + 20.272458*I
exp(p_1)-p_1 | %693 = 1.6045424 E-219 + 0.E-220*I \\ error is machine- epsilon
2-periodic points
That procedere for finding $1$-periodic points is of course well known. For the finding of $1$-periodic (or: fixed-) points we have also the LambertW-function and their branching behaviour for even an immediate solution. But none such function is known for the $n$-periodic points for $n\ge 2$. Thus we have to find our own procedure now...
For that search of $2$-periodic fixpoints I'll not go to explain the formulae in detail (it is trivial but too tedious for this answer-box), but I simply state that for the search for a solution to
$$ \exp(\exp(p_2) ) = p_2 \ne \exp(p_2) $$
it suffices to use -instead of the iterated exponential-function- the iterated branched logarithm-function, where the branch-index is allowed/required to change between the two steps of iteration:
$$ p_2 = \text{lb}( \text{lb}(p_2,k_0),k_1) $$
Initialize arbitrarily, iterate and, if desired, refine using the Newton-iteration on the $\exp(\exp())$ function.
Example:
p_2 = 1+I | %694 = 1 + I
p_2=lb(lb(p_2,-1),0) | %695 = 1.7063287 - 1.5078409*I
p_2=lb(lb(p_2,-1),0) | %696 = 1.9537456 - 1.4538908*I
p_2=lb(lb(p_2,-1),0) | %697 = 1.9430376 - 1.4429267*I
p_2= newton(p_2) | %698 = 1.9428329 - 1.4437804*I
exp(exp(p_2))-p_2 | %699 = 6.1238446 E-202 - 1.6330252 E-201*I
exp(p_2)-p_2 | %700 = -1.0588348 - 5.4784957*I
We have now:
p_2 | %701 = 1.9428329 - 1.4437804*I
exp(p_2) | %702 = 0.88399815 - 6.9222761*I
exp(exp(p_2)) | %703 = 1.9428329 - 1.4437804*I
which is indeed a $2$-periodic point.
It is meaningful now, to simplify the notation for finding $2$-periodic points. Just let us use a vectorial notation for the two branch-indexes and write
$$ p_2 = \text{Find}([k_1,k_2]) $$
implemented as Pari/GP-function (Pseudocode)
Find(K) = my(z=1+I,k_1=K[1],k_2=K[2]); \\ K is the vector [k_1,k_2]
for(t=1,5, \\ 5 or even only 3 iterations suffice to start Newton
z=lb(z,k_1);
z=lb(z,k_2);
);
return(Newton(z)) ;
Update A better routine, which also employs the Newton-iteration on the branched iterated logarithm (instead on iterated exponentiation):
\\Pari/GP
default(realprecision,200) \\ my usual numerical precision
pi2i=2*Pi*I \\ constant
{Find(K,maxerr=1e-100,maxit=25,z0=1+I)=my(err,L,z,l_prod,n=#K);
L=vector(n); \\ shall contain the sequence of logarithms/periodic points
\\ the following preconditioning is likely not needed at all
L[1]=z0; for(it=1,3,for(i=1,n,L[(i % n) + 1]=lb(L[i],K[i])));
z0=L[1];
\\ Newton-iteration on branched iterated logarithm,initial value z0
for(it=1,maxit,
z=z0;
L[1] = l_prod = z;
for(i=1,n, L[(i % n)+1]=z=log(z)+K[i]*pi2i; if(i<n,l_prod*=z) );
err = (z-z0)/(1/l_prod-1); \\ denominator contains derivative
z0 -= err;
if(abs(err) < maxerr,break());
);
return(L);}
end update
Using $Find([k_1,k_2])$ for $k_1=-3..3$ and $k_2=0..12$ I get the following chart of $1$- and $2$-periodic points (they are $1$-periodic when $k_1=k_2$) where the first parameter $k_1$ controls the color:

To see the $2$-periodicities, below is the same picture whith the pairs of periodic points connected by straight lines:
It is worth to note, that $2$-periodic points of the form $\text{Find}([k,-k])$ give pairs of complex conjugate numbers (big red diamonds), while $\text{Find}([k,k]) (=\text{Find}([k]))$ give $1$-periodic points (big brown circles). Note moreover that the $1$-periodic and the conjugate $2$-periodic points lay asymptotically on an (exponential) curve which can be seen, when the scale of the imaginary axis is taken as logarithmic (or even better transformed to $\sinh^{-1}()$ to see also the numbers with negative imaginary component ).
3-,4-,5-,... n-periodic points
If we extend the $\text{Find}([...])$ function to $3,4,5,...n$ entries $\in \mathbb Z$ in the vector-argument we find easily any $n$-periodic point which we like. $\text{Find}([0,0,1])$,$\text{Find}([0,0,-1])$,$\text{Find}([0,0,2])$,...$\text{Find}([0,1,1])$ ...
Finally $\text{Find}([k_1,k_2,k_3])$ with $k_1,k_2,k_3 \in \mathbb Z$ give all $3$-periodic points, and in case $k_1=k_2=k_3$ the $1$-periodic fixed points and in case $k_1=k_2=k_3=0$ the primary fixed point of the $\log()$-function. Unfortunately, as Y. Galidakis pointed out in his answer, the Newton-iteration for higher iterates of the $\exp()$ gets more&more involved; numerical checks up to $n=63$ and $20$ iterations in the $\text{Find}([...])$-function however looked very promising but are still under consideration.
See here some picture for examples for $n=3$,$n=5$,$n=11$,$n=31$ - periodic points.
Here the coloring is chosen to make the exemplars of one shape-family better discernible where one shape-family is meant as having chosen $k_1,k_2$ constant and only $k_3$ varying.

I have shown only one exemplar because the overlay of more exemplars of the shape-family makes the picture too chaotic
Here I show 4 exemplars of a very special shape-family by keeping all but one of the vectorial arguments to zero:
$$\text{Find}([0,0,...,0,k_{31}])$$
which produces the natural iteration-map for the $\log()$ for $30$ steps and then in one step adds $k_{31} \cdot C$ .
Of course, this typical shape-family exists analoguously for all $n$-periodic points.

Remark
This scheme, if it is really exhausting (what seems obvious to me), gives nice intuition in more general statements about the nature of the set-of-periodic points of the exponential function.
The number of $2$-periodic fixpoints is then that of $\mathbb Z^2$. All of them can be indexed like the rational numbers by a pair of 2 integer indexes.
The number of $n$-periodic fixpoints is then that of $\mathbb Z^n$. Indexing as before, but with $n$-tuple of integer indexes.
If I understand the topological concept of "dense subset" correctly it is immediately obvious, that the set of $n$-periodic points is "dense" (Shen/Rempe-Gillen give a reference to a proof)
Because there is so far no reason to assume that some class of $n$-periodic points might be missing/impossible, this answers one doubt in Galidakis' answer, by here claiming that $n$-periodic point for all $n$ exist and there are $\mathbb Z^n$ of them.(See also Shen/Rempe-Gillen, Theorem 1.1)
I don't know yet, whether it shall be possible to find, for instance for the $2$-periodic points, a simpler analytical description than this of the iterated branched logarithms, for instance by a parametrical definition of the curves on which that points lay, for instance to prove the exhaustiveness of my method for the set of 2-periodic orbits/points.
Shen, Zhaiming; Rempe-Gillen, Lasse, The exponential map is chaotic: an invitation to transcendental dynamics, Am. Math. Mon. 122, No. 10, 919-940 (2015). ZBL1361.37002.
Update: an article which deals with the question of $p_1$ (fixed-) points on the braches of the $\log()$-function is by Stanislav Sykora (2016) at his web-space here . Don't know really whether the exposition can be used at least as proof for the questions whether the set of fixed points $p_1$ found by this method here is really exhaustive, though.