Assume $n\ge m$, what is the estimation of
$$\sum_{k_1+\dots +k_m\,=\,n,\\ k_1\ge 1,\,\dots,\,k_m\ge 1} C_n^{k_1,\dots,k_m} \left(\frac{1}{k_1}+\frac{1}{k_2}+\dots +\frac{1}{k_m} \right)$$
where $C_n^{k_1,\dots,k_m}$ is the multinomial coefficient.
I am particularly interested in an upper bound, the tighter the better.
$\newcommand{\Z}{\mathbb Z}\renewcommand{\b}{\binom}$Here we will express the sum in question as an ordinary integral of an ordinary sum, which is much easier to analyze than the original expression -- which was a sum over an $m$-dimensional integer simplex, involving multinomial coefficients.
We will also analyze the simpler expression in the case when $m\to\infty$ and $n\ge cm\ln m$ for a real constant $c>1$. Then the resulting upper bound is optimal up to a constant factor.
For positive integers $n$ and $m$, let \begin{equation*} K^0:=\{k=(k_1,\dots,k_m)\in\Z_+^m\colon k_1+\dots+k_m=n\}, \end{equation*} where $\Z_+:=\{0,1,\dots\}$, and \begin{equation*} K:=\{k\in K^0\colon k_1\ge1,\dots,k_m\ge1\}. \end{equation*} Then, by the exchangeability of the $k_i$'s, the sum in question is \begin{equation*} S_{n,m}:=m\sum_{k\in K}\b nk\frac1{k_m}, \tag{0}\label{0} \end{equation*} where \begin{equation*} \b nk:=\b n{k_1,\dots,k_m}=\frac{n!}{k_1!\dots,k_m!}. \end{equation*} So, \begin{equation*} S_{n,m}=m\int_0^1\frac{dx}x\,S(K), \tag{1}\label{1} \end{equation*} where \begin{equation*} S(J):=\sum_{k\in J}\b nk x^{k_m} \end{equation*} for all $J\in P([m])$, where in turn $P([m])$ is the power set of the set $[m]:=\{1,\dots,m\}$. For $J\in P([m])$, let also \begin{equation*} K_J:=\{k\in K^0\colon k_j=0\ \forall j\in J\}, \end{equation*} so that $K_\emptyset=K^0$ and $K_{[m]}=\emptyset$ (since $n\ge1$). So, by the inclusion-exclusion principle, \begin{equation*} S(K)=\sum_{j=0}^{m-1}(-1)^j\sum_{|J|=j}S(K_J), \tag{2}\label{2} \end{equation*} where $\sum_{|J|=j}$ denotes the summation over all $J\in P([m])$ of cardinality $|J|=j$.
Next, for each $J\in P([m])$ with $|J|=j$ we have the following:
If $m\in J$, then $S(K_J)=(m-j)^n$.
If $m\notin J$, then $S(K_J)=(x+m-j-1)^n$.
Also, there are $\b{m-1}{j-1}$ sets $J\in P([m])$ with $|J|=j$ such that $m\in J$; and there are $\b{m-1}{j}$ sets $J\in P([m])$ with $|J|=j$ such that $m\notin J$. So, by \eqref{2}, \begin{align*} S(K)&=\sum_{j=0}^{m-1}(-1)^j\Big[\b{m-1}{j-1}(m-j)^n+\b{m-1}{j}(x+m-j-1)^n\Big] \\ &=\sum_{j=0}^{m-1}(-1)^j\b{m-1}{j}\big[(x+m-j-1)^n-(m-j-1)^n\big]; \end{align*} here we took into account that $\b{m-1}{-1}=0$ and $m-m=0$ and $n\ge1$. So, in view of \eqref{1}, \begin{equation*} S_{n,m}=m\sum_{j=0}^{m-1}(-1)^j\b{m-1}{j}I_{n,m}, \tag{3}\label{3} \end{equation*} where \begin{align*} I_{n,m}&:=\int_0^1\frac{dx}x\,\big[(x+m-j-1)^n-(m-j-1)^n\big] \\ &:=\int_0^1 dx\,\ln\frac1x\;n(x+m-j-1)^{n-1}, \end{align*} via integration by parts. So, by \eqref{3}, \begin{equation*} S_{n,m}=m\int_0^1 dx\,\ln\frac1x\;\sum_{j=0}^{m-1}(-1)^j\b{m-1}{j}n(x+m-j-1)^{n-1}. \tag{4}\label{4} \end{equation*}
The latter expression for $S_{n,m}$, which is an ordinary integral of an ordinary sum, is much easier to analyze than the original expression \eqref{0} for $S_{n,m}$ (which was a sum over an $m$-dimensional integer simplex, involving multinomial coefficients).
Specifics of analysis of the expression for $S_{n,m}$ in \eqref{4} depend on how $m$ varies with $n$.
For instance, suppose that
\begin{equation*}
n\ge m\ln m, \tag{5}\label{5}
\end{equation*}
which seems a natural assumption.
Then for
\begin{equation*}
a_j:=a_{n,m;j}(x):=\b{m-1}{j}n(x+m-j-1)^{n-1}
\end{equation*}
and $x\in(0,1)$ and $j\in\{0,\dots,m-2\}$ we have
\begin{equation*}
r_j:=\frac{a_{j+1}}{a_j}=\Big(1-\frac1{x+m-j-1}\Big)^{n-1}\frac{m-j-1}{j+1} \\
\le\Big(1-\frac1m\Big)^{n-1}(m-1)\le m e^{-n/m}\le1.
\end{equation*}
So, the $a_j$'s are nonincreasing in $j\in\{0,\dots,m-1\}$.
So, by \eqref{4},
\begin{align*}
S_{n,m}&=m\int_0^1 dx\,\ln\frac1x\;\sum_{j=0}^{m-1}(-1)^j a_{n,m;j}(x) \\
&\le m\int_0^1 dx\,\ln\frac1x\;a_{n,m;0}(x) \\
&=m\int_0^1 dx\,\ln\frac1x\;n(x+m-1)^{n-1}.
\tag{6}\label{6}
\end{align*}
To upper-bound the latter integral, take any $b\in(0,1)$.
Since $\ln\frac1x$ is decreasing and $n(x+m-1)^{n-1}$ is increasing in $x\in(0,1)$, by the Chebyshev integral inequality we have
\begin{align*}
&\int_0^b dx\,\ln\frac1x\;n(x+m-1)^{n-1} \\
& \le\frac1b\,\int_0^b dx\,\ln\frac1x\;\int_0^b dx\,n(x+m-1)^{n-1} \\
& =\ln\frac1b\;[(b+m-1)^{n}-(m-1)^{n}] \\
& \le\ln\frac1b\;(b+m-1)^{n} \\
& \le\ln\frac1b\;m^n\Big(1+\frac{b-1}m\Big)^{n} \\
&\le\ln\frac1b\;m^n e^{-(1-b)n/m}.
\tag{7}\label{7}
\end{align*}
Since $\frac{\ln(1/x)}{1-x}$ is decreasing in $x\in(0,1)$, with
\begin{equation*}
C_b:=\frac{\ln(1/b)}{1-b}
\end{equation*}
we get
\begin{align*}
&\int_b^1 dx\,\ln\frac1x\;n(x+m-1)^{n-1} \\
& \le C_b\int_b^1 dx\,(1-x)\;n(x+m-1)^{n-1} \\
& \le C_b\int_0^1 dx\,(1-x)\;n(x+m-1)^{n-1} \\
& = C_b\frac mn[m^n-(m-1)^n-n(m-1)^{n-1}]\le C_b\frac{m^{n+1}}n.
\tag{8}\label{8}
\end{align*}
Collecting \eqref{6}, \eqref{7}, and \eqref{8}, we get
\begin{align*}
S_{n,m}&\le B_{n,m}:=m^{n+1}\Big(\frac{C_b}n+\ln\frac1b\;e^{-(1-b)n/m}\Big).
\tag{9}\label{9}
\end{align*}
If we now strengthen condition \eqref{5} to \begin{equation*} n\ge cm\ln m \tag{5a}\label{5a} \end{equation*} for any real constant $c>\frac1{1-b}$, then for the bound $B_{n,m}$ in \eqref{9} we will have \begin{equation*} B_{n,m}\sim C_b \frac{m^{n+1}}n \end{equation*} as $m\to\infty$.
This may be compared with the following trivial bound on $S_{n,m}$: in view of \eqref{0}, \begin{equation*} S_{n,m}\le m\sum_{k\in K^0}\b nk=m^{n+1}. \tag{10}\label{10} \end{equation*}
So, given \eqref{5a}, the bound $B_{n,m}$ in \eqref{9} improves the trivial bound in \eqref{10} by a factor $\asymp n$. Following the lines of the above reasoning and noting that $\ln\frac1x>1-x$ for $x\in(0,1)$, one can see that this factor cannot be further improved, at least given \eqref{5a}.
Note also that, again assuming \eqref{5a}, the bound $B_{n,m}$ in \eqref{9} is better than the bound conjectured at the end of Brendan McKay's answer by a factor $\asymp m$.
Given that an upper bound is requested, here is one. By symmetry one can replace $\frac1{k_1}+\cdots+\frac1{k_m}$ by $\frac1{k_1}$ and multiply everything by $m$. Now do the sum over $k_1$ (which I'll call $k$) separately. $$\begin{align*}X &= m \sum_{k\ge 1} \frac 1k \sum_{k_2+\cdots+k_m=n;\, k_2\ge 1,\ldots,k_m\ge 1} \binom{n}{k,k_2,\ldots,k_m}\\ &= m \sum_{k\ge 1} \frac 1k \binom nk \sum_{k_2+\cdots+k_m=n;\, k_2\ge 1,\ldots,k_m\ge 1} \binom{n-k}{k_2,\ldots,k_m}. \end{align*}$$ The last sum has no closed form, but it is bounded above by dropping the conditions $k_2\ge 1,\ldots,k_m\ge 1$. Using the multinomial theorem, $$ X \le m \sum_{k\ge 1} \frac 1k \binom nk (m-1)^{n-k}.$$ This will be very accurate if $m$ is not too large, as the additional terms we added will be negligible then. Unfortunately we now have a hypergeometric function and I don't know of a sharp upper bound.
When $m$ is not very large the sequence is approximately gaussian with mode around $k=\frac nm$, so expanding around that point will provide an asymptotic expression. Continuing by guesswork, the sum should be about $\frac mn$ times the sum without the $\frac 1k$, which gives an estimate $$ \frac {m^2}{n} \bigl( m^n - (m-1)^n\bigr).$$ All this needs checking.
