This paper describes the development of an iterative three-dimensional parabolic equation solver that takes into account the effects of irregular boundaries and refraction from a layered atmosphere. A terrain-following coordinate transformation, based on the well-known Beilis-Tappert mapping, is applied to the narrow-angle parabolic equation in an inhomogeneous media. The main advantage of this approach, which has been used in two dimensions in the past, is the simplification of the impedance boundary condition at the earth's surface. The transformed initial-boundary value problem is discretized using the Crank-Nicholson marching scheme in the propagating direction and second-order finite-differences in the transversal plane. The proposed method relies on an efficient iterative fixed-point solver, which involves the inversion of tridiagonal matrices only. The accuracy of the method is evaluated through a comparison with boundary element simulations in a homogeneous atmosphere above a Gaussian hill. Results show that transversal scattering occurs in the shadow zone of the obstacle where the two-dimensional parabolic equation underestimates the pressure amplitude. The model is particularly suited for the simulation of infrasound in a three-dimensional environment with realistic topographies.