Context. Models of photon-dominated regions (PDRs) still fail to fully reproduce some of the observed properties, in particular the combination of the intensities of different PDR cooling lines together with the chemical stratification, as observed e.g. for the Orion Bar PDR. Aims. We aim to construct a numerical PDR model, KOSMA-\tau 3D, to simulate full spectral cubes of line emission from arbitrary PDRs in three dimensions (3D). The model is to reproduce the intensity of the main cooling lines from the Orion Bar PDR and the observed layered structure of the different transitions. Methods. We build up a 3D compound, made of voxels ("3D pixels") that contain a discrete mass distribution of spherical "clumpy" structures, approximating the fractal ISM. To analyse each individual clump the new code is combined with the KOSMA-\tau PDR model. Probabilistic algorithms are used to calculate the local FUV flux for each voxel as well as the voxel-averaged line emissivities and optical depths, based on the properties of the individual clumps. Finally, the computation of the radiative transfer through the compound provides full spectral cubes. To test the new model we try to simulate the structure of the Orion Bar PDR and compare the results to observations from HIFI/Herschel and from the Caltech Submillimetre Observatory (CSO). In this context new Herschel data from the HEXOS guaranteed-time key program is presented. Results. Our model is able to reproduce the line integrated intensities within a factor 2.5 and the observed stratification pattern within 0.016 pc for the [Cii] 158 \mu m and different 12/13 CO and HCO+ transitions, based on the representation of the Orion Bar PDR by a clumpy edge-on cavity wall. In the cavity wall, a large fraction of the total mass needs to be contained in clumps. The mass of the interclump medium is constrained by the FUV penetration. Furthermore, ...