A Fortran 77 subroutine is provided for implementing the Gibbs sampling procedure to a multiunidimensional IRT model for binary item response data with the choice of uniform and normal prior distributions for item parameters. In addition to posterior estimates of the model parameters and their Monte Carlo standard errors, the algorithm also estimates the correlations between distinct latent traits. The subroutine requires the user to have access to the IMSL library. The source code is available at http://www.siuc.edu/~epse1/sheng/Fortran/MUIRT/GSMU2.FOR. An executable file is also provided for download at http://www.siuc.edu/~epse1/sheng/Fortran/MUIRT/EXAMPLE.zip to demonstrate the implementation of the algorithm on simulated data.