The first-principle understanding of the processes in the Scrape-Off-Layer (SOL) of a tokamak is crucial for the developement of a thermonuclear reactor. Since the plasma temperature in the SOL is rather low, the plasma is typically not fully ionized, and the neutral atoms play an important role in determining the SOL regimes. We have derived a kinetic model for neutral atoms in the SOL that contains the fundamental elements of neutral dynamics, while remaining relatively simple. The model has been implemented in GBS [1], a three-dimensional numerical code developed to simulate SOL turbulence. GBS is based on the drift-reduced Braginiskii equations and is able to study the self-consistent formation of the plasma profile as the interplay of the plasma outflowing from the core, the parallel losses, and turbulent transport. Details of the numerical implementation are given and the first GBS results of coupled plasma turbulence and neutral dynamics are presented. [1] P. Ricci, et al., Plasma Phys. Control. Fusion 54 (2012) 124047