You realise foam has a built-in particle tracking layer? IMO the best thing to do would be to track the particles on the foam mesh, which automatically keeps track of cell locations and boundary collisions. If you simply have to calculate the particle trajectories with an external app, then you can map this back onto the foam particles before updating your tracking. You can also use the "interpolate" functions to do second order mapping of field values onto the particles.

I think the particle tracking layer wouldn't help, hence the coupled solvers should exchange their fluxes.

Selecting the cells by hand seems complicated for complex geometries. Maybe there is the possibility to define something like a region in the blockmeshdict-file, so that I can select the interesting area at once?