Theoretical atmospheric neutrino flux estimates serve as a crucial input for the determination of the neutrino mass hierarchy, the unitarity of the PMNS matrix and the atmospheric mixing angle $\theta_{23}$ in underground neutrino detectors, such as the Super-Kamiokande, IceCube DeepCore and KM3Net ORCA. With the expected reduction of detector-induced systematic uncertainties by the IceCube Upgrade, and the substantial gain in effective volume of the upcoming Hyper-Kamiokande and KM3NeT ORCA detectors, the theoretical uncertainty of the non-oscillated neutrino flux and flavor composition will ultimately impact the achievable precision of future measurements. In this work, we tackle the uncertainty associated with modeling of hadronic interactions, which has the largest effect on the calculation. We develop an empirical, data-driven model (DDM), derived from high-precision accelerator data from the recent CERN North Area (NA) fixed-target experiments, and a few simple model-dependent arguments. The model is well constrained in the intermediate energy range above a few GeV up to a hundred GeV and achieves good agreement with atmospheric muon data without explicitly using it. We compare our result to reference calculations of the atmospheric neutrino flux.