The progress in a 3-dimensional, non-local neoclassical transport simulation code "FORTEC-3D" is described. The main purpose of the code is to solve the drift-kinetic equation in general a 3-dimensional configuration using the δ f Monte Carlo method, and to calculate neoclassical fluxes and the time evolution of the ambipolar radial electric field simultaneously. This article explains new numerical schemes adopted in FORTEC-3D in order to overcome numerical problems, which happen especially in the cases where the bifurcation of radial electric field occurs. Examples of test simulation for an LHD magnetic field configuration with a bifurcated electric field are also shown. With improved numerical schemes, FORTEC-3D can calculate neoclassical fluxes and trace the time evolution stably for several ion collision times, which is sufficiently long to observe GAM damping and formation of the ambipolar electric field.