No, there is not a built-in function call for multivariate normal CDFs when d>2. For d=2 (bivariate), you can use the PROBBNRM function and for higher-dimensions you can use Monte Carlo estimation as follows:
/* Monte Carlo estimation for 3-D CDF of MVN with
CS covariance structure */
proc iml;
call streaminit(123);
d = 3; /* dimension */
N = 1e7; /* sample size */
mean = j(1, d, 0); /* mean of distribution */
Sigma = {1 0.5 0.5,
0.5 1 0.5,
0.5 0.5 1};
Z = randnormal(N, mean, I(d)); /* Z ~ MVN (0, I) */
v0 = {0 0 0}; /* cutoff values in each component (c,c,c) */
ComponentsInRegion = (Z < v0)[,+]; /* number of components in region */
group = (ComponentsInRegion=d); /* all 3 components in the region? */
MCEst = mean(group); /* proportion of obs in region */
print MCEst;