Probabilistic harmonic power flow (PHPF) in distribution networks with high penetration of photovoltaic generation and electric vehicle charging requires extensive Monte Carlo simulation, particularly when compliance must be evaluated under specific operating conditions. Conventional Monte Carlo analysis provides marginal output distributions, but these cannot be directly conditioned on individual operating regimes. We develop a conditional normalizing flow (CNF) surrogate that learns the full conditional distribution of per-bus harmonic voltages and total harmonic distortion indices given uncertain system states. This enables direct estimation of scenario-dependent compliance risk without rerunning simulations. The model is trained on 500,000 Monte Carlo harmonic simulations of the IEEE 33-bus distribution feeder using a conditional neural spline flow architecture. The surrogate achieves mean quantile errors of 0.005\% and 0.011\% at the 95th and 99th percentiles, respectively, a mean conditional exceedance error of 0.0009, and a probability integral transform deviation of 0.047, indicating accurate tail calibration. The results reveal strong regime dependence that is not visible in marginal statistics: at the most critical bus, the probability of exceeding the IEEE 519 limit of 8\% ranges from 0.095 under low-load conditions to 0.909 under high-load conditions, whereas the marginal estimate is 0.509. Once trained, the surrogate evaluates such conditional risk metrics in milliseconds through sampling, supporting rapid multi-scenario harmonic compliance analysis. This manuscript has been submitted to Electric Power Systems Research.