IQPACK: FORTRAN subroutines for the weights of interpolatory quadratures We present FORTRAN subroutines that implement the method described by the authors [Numer. Math. 40, 407-422 (1982; Zbl 0487.65014)] for the stable evaluation of the weights of interpolatory quadratures with prescribed simple or multiple knots. Given a set of knots and their multiplicities, the package generates the weights by using the zeroth moment $musb 0$ of w, the weight function in the integrand, and the (symmetric tridiagonal) Jacobi matrix J associated with the polynomials orthogonal on (a,b) with respect to w. There are utility routines that generate $musb 0$ and J for classical weight functions, but quadratures can be generated for any $musb 0$ and J supplied by the user. Utility routines are also provided that (1) evaluate a computed quadrature, applied to a user-supplied integrand, (2) check the polynomial order of precision of a quadrature formula, and (3) compute the knots and weights of simple Gaussian quadrature formula.