I recently had the same question. After emailing Gaussian support I was informed that the units are indeed $\sqrt{\text{km}/\text{mol}}$. From their reply I also believe that the derivatives are calculated with respect to the 3N-6 vibrational modes obtained through the application of the Sayvetz conditions as detailed in their whitepaper (Paragraph after Eqn. 11).
To understand why Gaussian uses these "odd" units, and how to get away from them, it's useful to look at the definition of the integrated absorption coefficient for the $i$th absorption band (peak) which corresponds to the $i$th vibrational mode
$$A_i = \frac{1}{nl}\int_{\text{band}}\ln\frac{I_0}{I}\text{d}\tilde{\nu} \ \ \ \ \ \ \ \text{[m mol}^{-1}] \tag{1}$$
which can be obtained directly from an experimental spectrum, and where $n$ is the concentration $[\text{mol m}^{-3}]$, $l$ the path length $[\text{m}]$, $I$ the transmitted intensity, $I_0$ the incident intensity, and $\tilde{\nu}$ wavenumber $[\text{m}^{-1}]$.
This can be calculated theoretically as
$$A_i = \frac{N_\mathrm{A}d_i}{12\epsilon_0 c^2 }\left|\frac{\partial\mathbf{\mu}}{\partial Q_i}\right|_0^2 \cdot \text{kg}(\text{u}) \cdot \text{D}^{2}\text{(C}^{2}\text{ Å}^{2}) \ \ \ \ \ \ \ \text{[m mol}^{-1}]\tag{2}$$
Where $N_\mathrm{A}$ is Avogradro's constant $[\text{mol}^{-1}]$, $d_i$ is the degeneracy of the $i$th mode, $\epsilon_0$ is the permittivity of free space $[\text{F m}^{-1}]$, $c$ is the speed of light $[\text{m s}^{-1}]$, while the "derivative term" contains a vector of the derivatives of the electric dipole with respect to the $x$, $y$, and $z$ components of the $i$th (mass-weighted) normal mode coordinate (evaluated at the equilibrium geometry) $[\text{D Å}^{-1}\text{u}^{-1/2}]$. I've included the final two terms as they explicitly detail the conversion factors required to arrive at the final units of $[\text{m mol}^{-1}]$
$$\text{kg}(\text{u}) = (1.661\times10^{-27})^{-1} \text{u kg}^{-1}\tag{3}$$
$$\text{D}^{2}\text{(C}^{2}\text{ Å}^{2}) = (3.336\times10^{-20})^{2} \text{ C}^{2} Å^{2} \text{ D}^{-2} \tag{4}$$
Evaluating the purely constant terms and conversion factors (precision truncated for brevity):
$$ \frac{N_\mathrm{A}}{12\epsilon_0c^2} \cdot \text{kg}(\text{u}) \cdot \text{D}^{2}\text{(C}^{2}\text{ Å}^{2}) \tag{5}$$
$$
\!\!\!\!\!\! = \frac{6.022\times10^{23}[\text{mol}^{-1}] \times (1.661\times10^{-27})^{-1} [\text{u kg}^{-1}] \times (3.336\times10^{-20})^{2} [\text{ C}^{2} Å^{2} \text{ D}^{-2}]}{12\times8.854\times10^{-12} [\text{F m}^{-1}]\times (2.998\times10^8 [\text{m s}^{-1}])^2} $$
$$
\!\!\!\!\!\!\!\!\!\!\!\!= 42.2561 \text{ km mol}^{-1} \text{u } Å^{2} \text{ D}^{-2} \tag{7}
$$
The Gaussian help pages give this exact same conversion factor, which (I believe) means that the dipole derivatives reported by Gaussian include the constant terms. To get the true derivatives in $[\text{D Å}^{-1}\text{u}^{-1/2}]$, you need to divide the reported derivatives in $[\text{km}^{1/2}\text{ mol}^{-1/2}]$ by $\sqrt{42.2561} \text{ km}^{1/2}\text{ mol}^{-1/2} \text{u }^{1/2} Å \text{ D}^{-1}$
This then means that the Intensity values printed by Gaussian $[\text{km mol}^{-1}]$ make a bit more sense, they're the sum of the squares of the dipole derivatives $[\text{km}^{1/2} \text{ mol}^{-1/2}]$ along $x$, $y$, and $z$.
Many textbooks and articles will state the $42.2561 \text{ km mol}^{-1} \text{u } Å^{2} \text{ D}^{-2}$ factor verbatim, and will often exclude the relevant conversion factors (even in explicit numerical evaluations), so be warned if you venture into the literature! As an example see Bernath's Spectra of Atoms and Molecules (3rd Ed.) Eqn. 7.246.
The above equations will need to be reworked if Beers Law is applied using the common logarithm rather than the natural logarithm, this amounts to a $1/\ln(10)$ term in the theoretical definition of $A_i$.
As has been said, you can get the derivatives without having to include that specific iop command, by looking in the formatted checkpoint (.fchk) file which can be obtained by running Gaussian's formchk utility on the checkpoint (.chk) file. However, these numbers don't match those in the .log file, as they are the derivatives with respect to single atom displacements along the x, y, and z axes. To obtain the derivatives with respect to displacements along the normal mode vectors, you first must generate the matrix $\mathbf{L}_{\text{CART}}$ detailed in the whitepaper - these are the displacement vectors printed by Gaussian (though in the .log file they are scaled by the square root of the reduced mass and so this scaling must be undone). Then for each mode, multiply each element of the displacement vector $\mathbf{L}_{\text{CART}}$ (corresponding to the motion of a single atom in x, y, or z in a given mode) by the the associated per-atom dipole derivative vector. Finally sum up all of the contributions to the derivative from each atom in each direction for a given mode and multiply by 31.2231 to convert from $[\text{e u}^{1/2}]$ to $[\text{ km}^{1/2}\text{ mol}^{-1/2}]$ and obtain the same values given by the IOp command.
