05AB1E, 26 23 22 bytes

-3 bytes thanks to Emigna

-1 byte thanks to Kevin Cruijssen

Lãœ.ΔIôDζ«D€í«ε€нÙgQ}P

Try it online!

R, 164 148 bytes

-many bytes thanks to Giuseppe.

n=scan()
`!`=function(x)sd(colSums(2^x))
m=function()matrix(sample(n,n^2,1),n)
while(T)T=!(l=m())|!(g=m())|!t(l)|!t(g)|1-all(1:n^2%in%(n*l+g-n))
l
g

Try it online!

Dramatically inefficient - I think it's even worse than other brute force approaches. Even for n=3, it will probably time out on TIO. Here is an alternate version (155 bytes) which works for n=3 in about 1 second.

Works by rejection. The function m draws a random matrix of integers between $1$ and $n$ (without forcing each integer to appear exactly $n$ times, or in different columns - this explains the slowness of the code, and is the only thing changed in the faster version). Call this function twice, to get the "latin" and "greek" squares l and g. Then we check:

1. 
   all(1:n^2%in%(n*l+g-n)): are there $n^2$ different
   pairs of values in l $times$ g?


2. are l and g latin squares?

To check that a matrix is a latin square, use the function !. Since point 1. above has been validated, we know that each integer appears $n$ times in each of l and g. Compute the column-wise sums of 2^l: these sums are all equal if and only if each integer appears once in each column (and the value of the sum is then $2^n+1-2$). If this is true for both l and its transpose t(l), then l is a latin square; same for g. The function sd, which computes the standard deviation, is an easy way to check whether all values of a vector are equal. Note that it doesn't work for the trivial cases $n=0$ and $n=1$ , which is OK according to the OP.

A final note: as often in R code golf, I used the variable T, which is initialized as TRUE, to gain a few bytes. But this means that when I needed the actual value TRUE in the definition of m (parameter replace in sample), I had to use 1 instead of T. Similarly, since I am redefining ! as a function different from negation, I had to use 1-all(...) instead of !all(...).

JavaScript (ES6),  159 147  140 bytes

Returns a flat array of $ntimes n$ pairs of non-negative integers.

This is a simple brute force search, and therefore very slow.

n=>(g=(m,j=0,X=n*n)=>j<n*n?!X--||m.some(([x,y],i)=>(X==x)+(Y==y)>(j/n^i/n&&j%n!=i%n),g(m,j,X),Y=X/n|0,X%=n)?o:g([...m,[X,Y]],j+1):o=m)(o=[])

Try it online! (with prettified output)

Commented

n => ( // n = input
 g = ( // g is the recursive search function taking:
 m, // m[] = flattened matrix
 j = 0, // j = current position in m[]
 X = n * n // X = counter used to compute the current pair
 ) => //
 j < n * n ? // if j is less than n²:
 !X-- || // abort right away if X is equal to 0; decrement X
 m.some(([x, y], i) => // for each pair [x, y] at position i in m[]:
 (X == x) + // yield 1 if X is equal to x OR Y is equal to y
 (Y == y) // yield 2 if both values are equal
 // or yield 0 otherwise
 > // test whether the above result is greater than:
 ( j / n ^ i / n && // - 1 if i and j are neither on the same row
 j % n != i % n // nor the same column
 ), // - 0 otherwise
 // initialization of some():
 g(m, j, X), // do a recursive call with all parameters unchanged
 Y = X / n | 0, // start with Y = floor(X / n)
 X %= n // and X = X % n
 ) ? // end of some(); if it's falsy (or X was equal to 0):
 o // just return o[]
 : // else:
 g( // do a recursive call:
 [...m, [X, Y]], // append [X, Y] to m[]
 j + 1 // increment j
 ) // end of recursive call
 : // else:
 o = m // success: update o[] to m[]
)(o = []) // initial call to g with m = o = []

Jelly,  21  20 bytes

-1 thanks to Nick Kennedy (flat output option allows a byte save of ż"þ`ẎẎQƑ$Ƈ $rightarrow$ F€p`Z€QƑƇ)

Œ!ṗ⁸Z€Q€ƑƇF€p`Z€QƑƇḢ

Try it online! (Too slow for 4 in 60s on TIO, but if we replace the Cartesian power, ṗ, with Combinations, œc, it will complete - although 5 certainly will not!)

How?

Œ!ṗ⁸Z€Q€ƑƇF€p`Z€QƑƇḢ - Link: integer, n
Œ! - all permutations of [1..n]
 ⁸ - chain's left argument, n
 ṗ - Cartesian power (that is, all ways to pick n of those permutations, with replacement, not ignoring order)
 Z€ - transpose each
 Ƈ - filter, keeping those for which:
 Ƒ - invariant under:
 Q€ - de-duplicate each
 F€ - flatten each 
 ` - use this as both arguments of:
 p - Cartesian product
 Z€ - transpose each
 Ƈ - filter, keeping those for which:
 Ƒ - invariant under: 
 Q - de-duplicate (i.e. contains all the possible pairs)
 Ḣ - head (just one of the Latin-Greaco squares we've found)

Haskell, 207 143 233 bytes

(p,q)!(a,b)=p/=a&&q/=b
e=filter
f n|l<-[1..n]=head$0#[(c,k)|c<-l,k<-l]$[]where
	((i,j)%p)m|j==n=[[]]|1>0=[q:r|q<-p,all(q!)[m!!a!!j|a<-[0..i-1]],r<-(i,j+1)%e(q!)p$m]
	(i#p)m|i==n=[[]]|1>0=[r:o|r<-(i,0)%p$m,o<-(i+1)#e(`notElem`r)p$r:m]

Try it online!

OK, I think I finally got it this time. It works fine for n=5, n=6 times out on TIO but I think that might just be because this new algorithm is INCREDIBLY inefficient and basically checks all possibilities until it finds one that works. I'm running n=6 on my laptop now to see if it terminates with some more time.

Thanks again to @someone for pointing out the bugs in my previous versions

C#, 520 506 494 484 bytes

class Pstatic void Main(string[]a)int n=int.Parse(a[0]);int[,,]m=new int[n,n,2];int i=n,j,k,p,I,J;R:for(;i-->0;)for(j=n;j-->0;)for(k=2;k-->0;)if((m[i,j,k]=(m[i,j,k]+ 1) % n)!=0)goto Q;Q:for(i=n;i-->0;)for(j=n;j-->0;)for(k=2;k-->0;)for(p=n;p-->0;)if(p!=i&&m[i,j,k]==m[p,j,k]for(i=n;i-->0;)for(j=n;j-->0;)System.Console.Write(m[i,j,0]+"-"+m[i,j,1]+" ");

The algorithm of findinf a square is very simple. It is... bruteforce. Yeah, it's stupid, but code golf is not about speed of a program, right?

The code before making it shorter:

using System;

public class Program

 static int[,,] Next(int[,,] m, int n)
 for (int i = 0; i < n; i++)
 
 for (int j = 0; j < n; j++)
 
 for (int k = 0; k < 2; k++)
 
 if ((m[i, j, k] = (m[i, j, k] + 1) % n) != 0)
 
 return m;
 
 
 
 
 return m;
 
 static bool Check(int[,,] m, int n)
 
 for (int i = 0; i < n; i++)
 
 for (int j = 0; j < n; j++)
 
 for (int k = 0; k < 2; k++)
 
 for (int p = 0; p < n; p++)
 
 if (p != i)
 if (m[i, j, k] == m[p, j, k])
 return false;
 
 for (int p = 0; p < n; p++)
 
 if (p != j)
 if (m[i, j, k] == m[i, p, k])
 return false;
 
 
 
 

 for (int i_1 = 0; i_1 < n; i_1++)
 
 for (int j_1 = 0; j_1 < n; j_1++)
 
 int i_2 = i_1;
 for (int j_2 = j_1 + 1; j_2 < n; j_2++)
 
 if (m[i_1, j_1, 0] == m[i_2, j_2, 0] && m[i_1, j_1, 1] == m[i_2, j_2, 1])
 return false;
 
 for (i_2 = i_1 + 1; i_2 < n; i_2++)
 
 for (int j_2 = 0; j_2 < n; j_2++)
 
 if (m[i_1, j_1, 0] == m[i_2, j_2, 0] && m[i_1, j_1, 1] == m[i_2, j_2, 1])
 return false;
 
 
 
 
 return true;
 
 public static void Main()
 
 int n = 3;
 Console.WriteLine(n);
 int maxi = (int)System.Math.Pow((double)n, (double)n*n*2);
 int[,,] m = new int[n, n, 2];
 Debug(m, n);
 do
 
 m = Next(m, n);
 if (m == null)
 
 Console.WriteLine("!");
 return;
 
 Console.WriteLine(maxi--);
 while (!Check(m, n));


 Debug(m, n);
 

 static void Debug(int[,,] m, int n)
 
 for (int i = 0; i < n; i++)
 
 for (int j = 0; j < n; j++)
 
 Console.Write(m[i, j, 0] + "-" + m[i, j, 1] + " ");
 
 Console.WriteLine();
 
 Console.WriteLine();
 

Now, if you want to test it with n=3 you wil have to wait like an hour, so here is another version:

public static void Main()

 int n = 3;
 Console.WriteLine(n);
 int maxi = (int)System.Math.Pow((double)n, (double)n*n*2); 
 int[,,] result = new int[n, n, 2];
 Parallel.For(0, n, (I) =>
 
 int[,,] m = new int[n, n, 2];
 for (int i = 0; i < n; i++)
 for (int j = 0; j < n; j++)
 
 m[i, j, 0] = I;
 m[i, j, 1] = I;
 
 while (true)
 
 m = Next(m, n);
 if (Equals(m, n, I + 1))
 
 break;
 
 if (Check(m, n))
 
 Debug(m, n);
 
 
 );

Update: forgot to remove "public".

Update: used "System." instead of "using System;"; Also, thanks to Kevin Cruijssen, used "a" instead of "args".

Update: thanks to gastropner and someone.

Wolfram Language (Mathematica), 123 bytes

P=Permutations
T=Transpose
g:=#&@@Select[T[Intersection[x=P[P@Range@#,#],T/@x]~Tuples~2,2<->4],DuplicateFreeQ[Join@@#]&]&

Try it online!

I use TwoWayRule notation Transpose[...,2<->4] to swap the 2nd and 4th dimensions of an array; otherwise this is fairly straightforward.

Ungolfed:

(* get all n-tuples of permutations *)
semiLSqs[n_] := Permutations@Range@n // Permutations[#, n] &;

(* Keep only the Latin squares *)
LSqs[n_] := semiLSqs[n] // Intersection[#, Transpose /@ #] &;

isGLSq[a_] := Join @@ a // DeleteDuplicates@# == # &;

(* Generate Graeco-Latin Squares from all pairs of Latin squares *)
GLSqs[n_] := 
 Tuples[LSqs[n], 2] // Transpose[#, 2 <-> 4] & // Select[isGLSq];

Python 3, 271 267 241 bytes

Brute-force approach: Generate all permutations of the pairs until a Graeco-Latin square is found. Too slow to generate anything larger than n=3 on TIO.

Thanks to alexz02 for golfing 26 bytes and to ceilingcat for golfing 4 bytes.

Try it online!

from itertools import*
def f(n):
 s=range(n);l=len
 for r in permutations(product(s,s)):
 if all([l(x[0]for x in r[i*n:-~i*n])*l(x[1]for x in r[i*n:-~i*n])*l(r[j*n+i][0]for j in s)*l(r[j*n+i][1]for j in s)==n**4for i in s]):return r

Explanation:

from itertools import * # We will be using itertools.permutations and itertools.product
def f(n): # Function taking the side length as a parameter
 s = range(n) # Generate all the numbers from 0 to n-1
 l = len # Shortcut to compute size of sets
 for r in permutations(product(s, s)): # Generate all permutations of all pairs (Cartesian product) of those numbers, for each permutation:
 if all([l(x[0] for x in r[i * n : (- ~ i) * n]) # If the first number is unique in row i ...
 * l(x[1] for x in r[i * n:(- ~ i) * n]) # ... and the second number is unique in row i ...
 * l(r[j * n + i][0] for j in s) # ... and the first number is unique in column i ...
 * l(r[j * n + i][1] for j in s) # ... and the second number is unique in column i ...
 == n ** 4 for i in s]): # ... in all columns i:
 return r # Return the square

Octave, 182 bytes

Brute force method, TIO keeps timing out and I had to run it a bunch of times to get output for n=3, but theoretically this should be fine. Instead of pairs like (1,2) it outputs a matrix of complex conjugates like 1+2i. This might be stretching the rule a bit, but in my opinion it stll fits the output requirements. There must be a better way to do the two lines under the functino declaration though, but I'm not sure at the moment.

function[c]=f(n)
c=[0,0]
while(numel(c)>length(unique(c))||range([imag(sum(c)),imag(sum(c.')),real(sum(c)),real(sum(c.'))])>0)
a=fix(rand(n,n)*n);b=fix(rand(n,n)*n);c=a+1i*b;
end
end

Try it online!