Let $A=[C_1,\cdots,C_n]\in M_n$ be a latin square, $M=A^TA$, $a_n=1\times n+2\times (n-1)+\cdots+n\times 1$, $b_n=\sum_{i=1}^ni^2=n(n+1)(2n+1)/6$ and $u=[1,\cdots,1]^T$. Note that $m_{i,j}=<C_i,C_j>\in[a_n,b_n]$ and $m_{i,i}=b_n$.
Moreover, for every $i$, $\sum_{j}m_{i,j}=n^2(n+1)^2/4$ and $c_n=\sum_{j|j\not= i}m_{i,j}=(n-1)n^2(n+1)(3n+2)/12$, which leads to the mean of the off-diagonal entries being $\mu_n=n(n+1)(3n+2)/12$.
Put $M=b_n.I+S$ where the diagonal of the symmetric matrix $S$ is $0$, $s_{i,j}\in [a_n,b_n]$ and $g_i(S)=\sum_{j|j\not= i}s_{i,j}-c_n=0$; we seek the max of $f:S\rightarrow \det(M)$ under the constraints $g_i(S)=0$. According to Lagrange, there is $(\lambda_i)_i\in\mathbb{R}^n$ s.t., for every $i\not= j$, $tr(adj(M)E_{i,j})+\lambda_i=0$, that is $z_{j,i}+\lambda_i=0=z_{i,j}+\lambda_i$ where $z_{i,j}$ is the $i,j$ entry of $adj(M)$; then, for every $i$, the $(z_{i,j})_{j|j\not= i}$ are equal and, since $adj(M)$ is symmetric, the $(z_{i,j})_{i\not= j}$ are equal. Finally, $M^{-1}$ is in the form $D+\alpha U$ where $D=diag(d_i)$ is diagonal, $\alpha\in\mathbb{R}$ and $U=[u_{i,j}]$ is defined by $u_{i,i}=0$ and , if $i\not= j$, then $u_{i,j}=1$. Therefore $MM^{-1}=(b_nI+S)(D+\alpha U)=b_nD+SD+\alpha b_n U+\alpha SU=I$. We consider the diagonals of RHS and LHS: for every $i$, $b_nd_i+\alpha c_n=1$; consequently, the $(d_i)_i$ are equal and $M^{-1}$ is in the form $M^{-1}=\beta I+\alpha U$. Thus $S$ is a polynomial in $U$ and all the non-zero entries of $S$ are equal to $\mu_n$.
Conclusion. The (ideal) maximum of $|\det(A)|$ is reached when the $(m_{i,j})_{i\not= j}$ are $\mu_n$. Note that, in general, $\mu_n$ is not an integer and the considered ideal maximum is not reached; yet, the idea is to find a matrix $A$ s.t. $A^TA$ is close to $b_nI+\mu_n U$.
The case $n=5$. Then $\mu_5=42.5$ and $\sqrt{\det(b_5I+\mu_5 U)}\approx 2343.75$. The real maximum $2325$ is reached (according to Peter) with a circulant matrix s.t. the $(m_{i,j})_{i\not= j}$ are $42$ or $43$.
The case $n=6$. Then $\mu_6=70$ but $\sqrt{\det(b_6I+\mu_6 U)}\approx 42439.23$ is not an integer. A candidate for the real maximum is $41895$; it is reached with a circulant matrix s.t. the $(m_{i,j})_{i\not= j}$ are $69,70,71$.
The case $n=9$. Then $\mu_9=217.5$ and $\sqrt{\det(b_9I+\mu_9 U)}\approx 934173632.8$. A candidate for the real maximum is $929587995$; it is reached with a circulant matrix and a sudoku; the last one is s.t. the $(m_{i,j})_{i\not= j}$ are $216,217,218,219$.
Remark. The above reasoning works, in the same way and with the same conclusion, when we replace $M$ with $M'=AA^T$. Thus, for such an ideal solution, $A^TA=AA^T$, that is, $A$ is a normal matrix. Thus, good candidates for the Graal are the circulant matrices -that are normal-.
EDIT. We consider the Cholesky decomposition $b_nI+\mu_n.U=C^TC$; then $A^TA=b_n+\mu_n.U$ iff $A=PC$ where $P\in O(n)$. Finally, the idea is to rotate the known triangular matrix $C$ to make it close to an integer matrix (of course, it is not easy). For example, when $n=9$, a candidate solution $B$ has $[1,2,4,3,8,9,5,6,7]$ as first row; clearly, $B=RC$, where $R$ is close to an orthogonal matrix $T$ (consider the polar decomposition of $R$); if we replace $R$ with $T$, then we obtain $[1.02,2.07,4.08,2.89,7.99,9.02,4.95,5.92,7.07]$ as new first row; of course, this is this last form that must be discovered. We may also require that $A$ is normal, that is equivalent to the linear equation in the unknown $P$: $MP-PCC^T=0$. Unfortunately, the associated kernel is big; for example,when $n=9$, its dimension is $65$.
Best Answer
This may provide a little more insight into your question. It shows if you have an n by n matrix and n is even, you can produce singular matrices which, by definition, all have a determinant of 0.
Here is a paper you may want to look at, published in the Mathematical Association of America's Journal (Amer. Math. Monthly, Nov. 2000):
Even Order Regular Magic Squares are Singular, R. Bruce Mattingly