A 1k-dimensional multivariate normal distribution is made discrete by partitioning the k-dimensional Euclidean space with rectangular grids. The collection of probability integrals over the partitioned cubes is a k-dimensional contingency table with ordered categories. It is shown that loglinear model with main effects plus two-way interactions provides an accurate approximation for the k-dimensional table. The complete multivariate normal integral table is computed via the iterative proportional fitting algorithm from bivariate normal integral tables. This approach imposes no restriction on the correlation matrix. Comparisons with other numerical integration algorithms are reported. The approximation suggests association models for discretized multivariate normal distributions and contingency tables with ordered categories.