First of all,
let's see what the result should be
for general $k$.
$\begin{array}\\
f(k, n)
&=\dfrac{k^{kn} (n!)^k}{(kn)!}\\
&\approx \dfrac{k^{kn} (n^n\sqrt{2\pi n}e^{-n})^k}{(kn)^{kn}\sqrt{2\pi kn}e^{-kn}(kn)!}\\
&= \dfrac{k^{kn} n^{kn}(2\pi)^{k/2} n^{k/2}e^{-kn}}{k^{kn}n^{kn}\sqrt{2\pi kn}e^{-kn}}\\
&= \dfrac{(2\pi)^{(k-1)/2} n^{(k-1)/2}}{\sqrt{ k}}\\
\end{array}
$
As a check,
this gives
$f(3, n)
\approx \dfrac{2\pi n}{\sqrt{ 3}}
$
which agrees with your result,
since
$(k-1)/2 = 1$ for $k=3$.
Note that the
$(3n+1)!$ was sort of a fake,
since
$3n(3n)!$
would have given the same limit.
If you want to so
some sort of
Riemann sum,
you would have to take logs
and use
$\ln(n!)
=\sum_{i=1}^n \ln(i)
\approx \int_{i=1}^n \ln(x)dx
=(x \ln(x)-x)\big|_0^n
=n\ln(n)-n
$.
This is quite close
to Stirling's
$\ln(n!)
\approx n\ln(n)-n+\frac12(\ln(n)+\ln(2\pi))
$.
The reason is that
if $f$ is monotonic
then the error in using
$\sum_{i=1}^n f(i)
\approx \int_{i=1}^n f(x)dx
$
is bounded by
$f(0)+f(n)$
which is, in this case,
$\ln(n)$.
This can be proved
by using
$\min(f(n), f(n+1))
\le \int_n^{n+1} f(x) dx
\le \max(f(n), f(n+1))
$.
What you get
for $\ln f(k, n)=g(k, n)$
is
$\begin{array}
fg(k, n)
&\approx kn\ln(k)+k\ln(n!)-\ln((kn)!)\\
&=kn\ln(k)+k(n\ln(n)-n)-(kn\ln(kn)-kn)\\
&=kn\ln(k)+kn\ln(n)-kn-(kn(\ln(k)+\ln(n))-kn)\\
&=kn\ln(k)+kn\ln(n)-kn-(kn(\ln(k)+\ln(n))-kn)\\
&=0\\
\end{array}
$
The reason for this
is that the error
in
$\ln(n!)
\approx n \ln(n)-n
$
is of order
$\ln(n)$,
so the error in
$\ln(n!^k)$
is of order
$k\ln(n)
=\ln(n^k)
$
which fits the precise result
nicely.
In other words,
you would have to use
a more precise approximation
to $\sum_{i=1}^n \ln(i)$
and the best you could do
would be the same as
Stirling's approximation
with the constant
($\sqrt{2\pi}$)
being undetermined.
More precisely,
if you used
$\ln(n!)
\approx n\ln(n)-n+\frac12 \ln(n)+c
$,
you would get
$\ln(f(k, n))
\approx (k-1)c+\frac{k-1}{2}\ln(n)-\frac12 \ln(k)
$.
And that's about all I can think of.