A simple way to analyze the adiabatic lapse rate is to consider removing the parcel from its environment, changing its height $z$ in a gravitational field (corresponding to a potential energy $mgz$, with mass $m$ and gravitational acceleration $g$), and reinserting it into the environment.
In the absence of heat transfer (i.e., adiabatic conditions), the energy required to create a system at temperature $T$ and insert it into an environment—pushing the surroundings at pressure $P$ out of the way—is the enthalpy $H(T,P)$. Any increase in the potential energy is "paid for" by this enthalpy.
From this, we find that $H(T,P)+mgz$ is a constant for such adiabatic motion, i.e., $dh=-g\,dz$, with specific enthalpy $h$.
The change in specific enthalpy of an ideal gas is simply $dh=c_PdT$, so the dry adiabatic lapse rate is $$\Gamma_\text{dry}\equiv\left(\frac{dT}{dz}\right)_\text{dry}=-\frac{g}{c_P}.$$
(Sometimes the lapse rate is defined as the negative of this value; in this case, switch the signs here and below.)
Propagation of error tells us that for small changes in $g$ or $c_P$, $\Gamma_\text{dry}$ changes by
$$\frac{\Delta\Gamma_\text{dry}}{\Gamma_\text{dry}}=\left(\frac{1}{g_0}\right)\Delta g-\left(\frac{1}{c_{P,0}}\right)\Delta c_P.$$
Let's now consider a combination of ideal gases $i$. The change in enthalpy is $dH=\sum_i m_ic_{P,i}dT$:
$$\Gamma_\text{dry, mixture}\equiv\left(\frac{dT}{dz}\right)_\text{dry, mixture}=-\frac{g\sum_im_i}{\sum_im_ic_{P,i}}.$$
How about a combination of ideal gases $i$ saturated with a combination of vapors $j$ and suspended condensed matter $k$? The enthalpy is $H=H_0+\sum_i m_ic_{P,i}T+\sum_j m_j(c_{P,j}T+L_j)+\sum_k m_kc_{P,k}T$, where $L_j$ is the latent heat for condensation of vapor $j$. (For simplicity, I'm not considering successive phase changes, e.g., melting and then freezing. You can extend as you like.)
To recover the equation you show, assume that minimal water condenses or that precipitation falling from the parcel only minimally affects the enthalpy (sometimes called the pseudoadiabatic assumption), that only air and water vapor exist (total mass $m$, vapor mass $m_\text{vapor}$), and that the humidity only minimally affects the air specific heat capacity $c_P$:
$$H=H_0+mc_{P}T+m_\text{vapor}L;\tag{enthalpy expression}$$
$$dH=mc_{P}dT+dm_\text{vapor}L;\tag{differentiating}$$
$$dh(=-g\,dz)=c_{P}dT+\frac{dx_\text{vapor}}{dT}L\,dT=\left(c_{P}+\frac{dx_\text{vapor}}{dT}L\right)dT;\tag{per mass}$$
$$\frac{dT}{dz}=-\frac{g}{c_{P}+\frac{dx_\text{vapor}}{dT}L};\tag{rearranging}$$
where $x_\text{vapor}\equiv\frac{m_\text{vapor}}{m}$, which scales with the vapor partial pressure as $\frac{p_\text{vapor}}{P}$, so
$$\frac{dx_\text{vapor}}{dT}=\frac{x_\text{vapor}}{p_\text{vapor}}\frac{dp_\text{vapor}}{dT}-\frac{x_\text{vapor}}{P}\frac{dP}{dT};\tag{differentiating}$$
$$\frac{dx_\text{vapor}}{dT}=\frac{x_\text{vapor}}{p_\text{vapor}}\frac{dp_\text{vapor}}{dT}-\frac{x_\text{vapor}}{P}\frac{dP}{dz}\frac{dz}{dT};\tag{expanding}$$
$$\frac{dx_\text{vapor}}{dT}=\frac{x_\text{vapor}L}{T^2R_\text{vapor}}-\frac{x_\text{vapor}}{P}\frac{dP}{dz}\frac{dz}{dT};\tag{Clausius–Clapeyron}$$
$$\frac{dx_\text{vapor}}{dT}=\frac{x_\text{vapor}L}{T^2R_\text{vapor}}+\frac{x_\text{vapor}g}{TR_\text{air}}\frac{dz}{dT},\tag{barometric}$$
with some gas constants $R$, which ultimately reduces to
$$\Gamma_\text{saturated}\equiv\left(\frac{dT}{dz}\right)_\text{saturated}=-g\frac{1+\frac{x_\text{vapor}L}{TR_\text{air}}}{c_P+\frac{x_\text{vapor}L^2}{T^2R_\text{vapor}}},$$
which matches the form of your equation.
