Computational Fluid Dynamics (CFD) model was developed for the simulation of the neutrally-stratified atmospheric boundary layer flow over hills, whereas the additional windsource term was implemented into the steady Reynolds-averaged Navier Stokes (RANS) momentum equation. Such an additional volume force ensured computational shear stress distribution in agreement with the theoretical values accompanied with the negligible pressure gradient along the computational domain. The wind-source term was additionally corrected during the precursor periodic two- dimensional simulation in an empty computational domain in order to reduce discrepancies between computational and theoretical shear stress distributions. The calibration of the wind-source term allowed for achieving a good agreement between the computational and theoretical mean velocity values, whereas a periodic precursor domain technique also ensured flow homogeneity. The mean velocity, turbulence kinetic energy and dissipation rate profiles were used as an inflow boundary condition for the simulations of the atmospheric boundary layer (ABL) flow over a generic two-dimensional cosine-shaped hill for three hill slopes. The obtained computational results indicate that the proposed model is able to provide satisfactory results for the flow over hilly terrains because the computed flow separation and reattachment points agree well with previously reported results.