The exact number of latin squares of order n is known only up to n = 11, and these counts are reflected in the OEIS and Wikipedia entries already linked in the Comments.
For counting "diagonal" latin squares I could not find any ready online reference, so I wrote a quick Prolog program that provides solutions for the smaller orders.
Note that for the sake of counting we may assume a particular first row of the latin square, so that if the first row is $(1,2,...,n)$ it is in partially reduced form (reduced form having also the first column in standard order).
For n = 1 there is the one (trivial) case.
For n = 2,3 one quickly finds there are no diagonal latin squares.
For n = 4 there are two cases in partially reduced form. Therefore the total count is 2*4! = 48.
For n = 5 there are eight cases in partially reduced form, for a total count of 8*5! = 960.
My code is not optimized for speed but is clear, so the next step is to refactor and impose constraints in a more efficient manner to make some higher orders feasible.
Added:
For n = 6 there are 128 cases in partially reduced form, for a total count of 128*6! = 92160. This result was confirmed by a computation using a different normalization, see discussion below.
For n = 7 my preliminary result is a total count of 171200*7! = 862,848,000 solutions.
Using the partially reduced form (first row in standard order) the n = 6 case was enumerated in 6.3 seconds. Analyzing the time consumed in building stages of the search space suggested that normalizing the main diagonal (in standard order) would be better. I posted a Question here, Counting some special derangements, motivated by how this constrains the choice of anti-diagonals. Another feature of this approach is cutting the enumeration in half via the action of matrix transpose.
These efficiencies reduced the computation time for n = 6 to 1.5 seconds. Despite their use in the n = 7 case, the increased size of the search space led to a running time of 1776 seconds, nearly half an hour.
Such a three-order of magnitude increase in running time suggests that the n = 8 case will not be easily solved on a single CPU, even if the code were compiled instead of interpreted. [Update: For the n = 8 case there are 7450228224*8! = 300,393,201,991,680 diagonal latin squares. This was computed on a single CPU through a combination of speeding up the code that searches on a single main diagonal/anti-diagonal assignment and reducing the number of such assignments by equivalences. The code now runs such a single assignment in 1-2 hours. Given normalization of the main diagonal to the identity sequence $1,2,..,8$, there are 4,752 possibilities (µ-derangements) for the anti-diagonal. A first reduction to 630 equivalences classes is made by using orbits under a dihedral group action. A further reduction to 114 enlarged equivalence classes results from taking "orbits of orbits" under a group isomorphic to $(\mathbb{Z}/2\mathbb{Z})^4$ acting by conjugation. Finally it appears we actually need only count 18 distinct cases due to equivalences of a more contingent nature than those group symmetries. Similar results should apply to solving n = 9, though searching each given µ-derangement will take more time.]
While seeking citations for counts of diagonal latin squares I came across this interesting 2009 paper by Zhang and Ma, Counting solutions for the N-queens and Latin-square problems by Monte Carlo simulations, which is freely available through NIH's Pub Med Central. Perhaps their technique could be adapted for estimating counts of diagonal latin squares of "large" orders. More recently (2017) a report Using Volunteer Computing to Study Some Features of Diagonal Latin Squares describes tactics to count them up to order $n=8$.
This question has a definite answer, which is to use the Polya Enumeration Theorem, which is more powerful than Burnside and includes it as a special case. It will produce a multivariate generating function for all distributions of colors, from which you may extract the desired coefficient for a given configuration. You compute the cycle index of the face permutation group $G$ of the polyhedron and evaluate it using the substitution
$$a_k = X_1^k + X_2^k + \cdots + X_n^k$$
if you have $n$ colors. Then you extract the coefficient for your particular configuration, which is written like this:
$$[X_1^2 X_2^2 \cdots X_n^2] Z(G)(X_1+X_2+\cdots+X_n).$$
Here are two simple examples where we restrict ourselves to rigid motions (no reflections). Suppose we have a tetrahedron whose faces we want to color. The cycle index of $G$ for this case contains the identity, which contributes $$a_1^4.$$ Then there are rotations about axes passing through the middle of a face and the opposing vertex (four such axes), which fix that face and contribute $$4\times 2\times a_1 a_3.$$ Finally there are three flips by 180 degrees about an axis passing through the midpoints of two vertex-disjoint edges (three such axes), which contribute $$3\times a_2^2.$$ This gives the cycle index
$$Z(G) = \frac{1}{12} (a_1^4 + 8 a_1 a_3 + 3 a_2^2).$$
We can recover Burnside from this via
$$Z(G)(X_1+X_2+\cdots+X_n)_{X_1=X_2=\cdots=X_n=1}$$
which gives the formula
$$\frac{1}{12} (n^4 + 8 n^2 + 3 n^2) = \frac{1}{12} n^4 + \frac{11}{12} n^2.$$
for colorings of the tetrahedron with at most $n$ colors.
This sequence is documented as OEIS A006008.
Now to answer your question, suppose we have three colors $R$, $G$ and $B$ and we want the number of colorings where two colors are used twice. This is where the substituted cycle index comes into play. Making the substitution as described above we have
$$Z(G)(R+B+G)\\ =
1/12\, \left( R+G+B \right) ^{4}+2/3\, \left( R+G+B \right) \left( {R}^{3
}+{G}^{3}+{B}^{3} \right) +1/4\, \left( {R}^{2}+{G}^{2}+{B}^{2} \right) ^{
2}$$
which expands to
$${B}^{4}+{B}^{3}G+{B}^{3}R+{B}^{2}{G}^{2}+{B}^{2}GR+{B}^{2}{R}^{2}+B{G}^{3}
+B{G}^{2}R\\+BG{R}^{2}+B{R}^{3}+{G}^{4}+{G}^{3}R+{G}^{2}{R}^{2}+G{R}^{3}+{R}
^{4}.$$
Extracting coefficients we see that
$$[R^2 G^2]Z(G)(R+G+B)+[R^2 B^2]Z(G)(R+G+B)+[G^2 B^2]Z(G)(R+G+B) = 3,$$
so there are three colorings up to rotation of the faces of the tetrahedron with three available colors using two colors twice.
As a second example consider painting the faces of a cube with three colors. We need another cycle index call it $Z(H)$, the one of the face permutation group of the cube, which is actually a fairly standard computation. There is the identity, which contributes $$a_1^6.$$ There are rotations about an axis passing through opposite faces (three such axes) which fix those faces, giving $$3\times (2 a_1^2 a_4 + a_1^2 a_2^2).$$ There are rotations about axes passing through opposite vertices (four such axes), giving $$4\times 2 \times a_3^2.$$ Finally there are rotations about axes passing through the midpoints of pairs of opposite edges (six such pairs), which contribute $$6\times a_2^3.$$
This gives the cycle index
$$Z(H) = \frac{1}{24}
(a_1^6 + 6 a_1^2 a_4 + 3 a_1^2 a_2^2 + 8 a_3^2 + 6 a_2^3).$$
We recover Burnside from this as above, getting
$$\frac{1}{24} (n^6 + 6 n^3 + 3 n^4 + 8 n^2 + 6 n^3)
= \frac{1}{24} (n^6 + 3 n^4 + 12 n^3 + 8 n^2)
\\= \frac{1}{24} n^6 + \frac{1}{8} n^4 + \frac{1}{2} n^3 + \frac{1}{3} n^2.$$
This is sequence is OEIS A047780.
Now to get the colorings using each color twice we compute the substituted cycle index
$$Z(H)(R+G+B)=
1/24\, \left( R+G+B \right) ^{6}+1/4\, \left( R+G+B \right) ^{2} \left( {R
}^{4}+{G}^{4}+{B}^{4} \right) \\+1/8\, \left( R+G+B \right) ^{2} \left( {R}^
{2}+{G}^{2}+{B}^{2} \right) ^{2}+1/3\, \left( {R}^{3}+{G}^{3}+{B}^{3}
\right) ^{2}+1/4\, \left( {R}^{2}+{G}^{2}+{B}^{2} \right) ^{3}$$
which expands to
$${B}^{6}+{B}^{5}G+{B}^{5}R+2\,{B}^{4}{G}^{2}+2\,{B}^{4}GR+2\,{B}^{4}{R}^{2}
+2\,{B}^{3}{G}^{3}\\+3\,{B}^{3}{G}^{2}R+3\,{B}^{3}G{R}^{2}+2\,{B}^{3}{R}^{3}
+2\,{B}^{2}{G}^{4}+3\,{B}^{2}{G}^{3}R+6\,{B}^{2}{G}^{2}{R}^{2}\\+3\,{B}^{2}G
{R}^{3}+2\,{B}^{2}{R}^{4}+B{G}^{5}+2\,B{G}^{4}R+3\,B{G}^{3}{R}^{2}+3\,B{G}
^{2}{R}^{3}+2\,BG{R}^{4}\\+B{R}^{5}+{G}^{6}+{G}^{5}R+2\,{G}^{4}{R}^{2}+2\,{G
}^{3}{R}^{3}+2\,{G}^{2}{R}^{4}+G{R}^{5}+{R}^{6}.$$
Extracting coefficients we now see that
$$[R^2 G^2 B^2]Z(H)(R+G+B) = 6,$$
so there are six colorings using each of the three colors exactly twice.
The following MSE link points to a chain of similar computations.
Best Answer
Let's consider separately cases where there are $4,2$ and $3$ different colours in the first column.
Case 1:
$4$ different colors. This has only one pattern since the outermost colors must switch places with the center two colors and there is only one way this can happen. Colors can be chosen in $4!=24$ different ways. $$\begin{bmatrix}a&c&a&c&a\\b&d&b&d&b\\c&a&c&a&c\\d&b&d&b&d\end{bmatrix}$$
Case 2:
$2$ different colors. First column can be constructed in $4\cdot3=12$ ways. There can only be $2$ colors per column, but for every column from $2$ to $5$ there are two possible arrangements. In total there are $12\cdot2^4=192$ solutions.
$$\begin{bmatrix}a&c&a&c&a\\b&d&b&d&b\\a&c&a&c&a\\b&d&b&d&b\end{bmatrix}$$
Case 3:
3 different colors.
I) The first column can not be of the form $$\begin{bmatrix}a\\b\\c\\a \end{bmatrix}$$ (this would imply the two center colors would both be $d$)
II) If the same colors are arranged as follows there is only one pattern (since third row is fixed):
$$\begin{bmatrix}a&d&a&d&a\\b&c&b&c&b\\a&d&a&d&a\\c&b&c&b&c\end{bmatrix}$$
The colors can be chosen in $4\cdot3\cdot2=24$ ways.
III) The last possibility is case II) flipped. $24$ more combinations.
So in total there are $$24+192+2\cdot 24=264$$ ways color the squares such that the conditions are fulfilled.