In this paper Fortran subroutines for the evaluation of the discrete form of the Helmholtz integral operators L_k, M_k, M_k^t and N_k for two-dimensional, three-dimensional and three-dimensional axisymmetric problems are described. The subroutines are useful in the solution of Helmholtz problems via boundary element and related methods. The subroutines have been designed to be easy to use, reliable and efficient. The subroutines are also flexible in that the quadrature rule is defined as a parameter and the library functions (such as the Hankel, exponential and square root functions) are called from external routines. The subroutines are demonstrated on test problems arising from the solution of the Neumann problem exterior to a closed boundary via the Burton and Miller equation.