def METHOD_NAME(z_above_reservoir, rhod):
    if z_above_reservoir < 0:
        return 0
    water_vapour_mixing_ratio = self.water_vapour_mixing_ratio(
        z_above_reservoir
    )
    d_water_vapour_mixing_ratio__dz = derivative(
        self.water_vapour_mixing_ratio, z_above_reservoir
    )
    T = self.formulae.state_variable_triplet.T(
        rhod[0], self.thd(z_above_reservoir)
    )
    p = self.formulae.state_variable_triplet.p(
        rhod[0], T, water_vapour_mixing_ratio
    )
    lv = self.formulae.latent_heat.lv(T)
    return self.formulae.hydrostatics.drho_dz(
        g, p, T, water_vapour_mixing_ratio, lv
    ) / (
        1 + water_vapour_mixing_ratio
    ) - rhod * d_water_vapour_mixing_ratio__dz / (
        1 + water_vapour_mixing_ratio
    )