Multiscale random Bernstein polynomial (msBP) priors exhibit many favorable properties for nonparametric Bayesian inference. In msBP, an infinite tree of probability weights is generated from a generalization of the stick-breaking process representation of the Dirichlet process, with each tree scale featuring a weighted random Bernstein polynomial. The prior can generate densities with locally varying smoothness, accommodating abrupt local changes. The degree of smoothness of the resulting function approximation is determined by a hyperparameter that controls the decline in probabilities over the scales. In this paper we extend the multiscale Bernstein polynomial prior from [0, 1] analyzed in the previous literature to the unit hypercube \([0,1]^d\) for the link function in copula dependence modeling, accounting for the uncertainty inherent in the degree of approximation smoothness by Bayesian Model Averaging over the smoothness hyperparameter. The favorable properties of the resulting non-parametric copula model are attained at the cost of an increased computational burden in practical implementation. We provide details of the implementation algorithm that is based on large-scale parallelization tailored to the internal model structure. We implement the algorithm on distributed nodes of a Unix High Performance Computing system with Graphical Processing Units (GPU) acceleration via a Message Passing Interface (MPI) in Modern Fortran. The algorithm scales well across MPI ranks. The routines can be pre-compiled into numerical libraries and invoked from high-level languages such as Python or R.