A new three-dimensional code, VENUS+f, for neoclassical transport calculations in nonaxisymmetric toroidal systems is presented. Numerical drift orbits from the original VENUS code and the f method developed for tokamak transport calculations are combined. The first results obtained with VENUS+f are compared with neoclassical theory for different collisional regimes in a JT-60 tokamak test case both for monoenergetic particles and for a Maxwellian distribution; good agreement is found. Successful benchmarking of the bootstrap current in the Wendelstein 7-X configuration with the DKES code for different collisionality regimes as well as further VENUS+f developments are described.