Recent advances in global ocean prediction systems are fostered by the needs of accurate representation of mesoscale processes. The day-by-day realistic representation of its variability is hampered by the scarcity of observations as well as the capability of assimilation systems to correct the ocean states at the same scale. This work extends a 3DVAR system designed for oceanic applications to cope with global eddy-resolving grid and dense observational datasets in a hybridly parallelized environment. The efficiency of the parallelization is assessed in terms of both scalability and accuracy. The scalability is favored by a weak-constrained formulation of the continuity requirement among the artificial boundaries implied by the domain decomposition. The formulation forces possible boundary discontinuities to be less than a prescribed error and minimizes the parallel communication relative to standard methods. In theory, the exact solution is recovered by decreasing the boundary error toward zero. In practice, it is shown that the accuracy increases until a lower bound arises, because of the presence of the mesh and the finite accuracy of the minimizer. A twin experiment has been set up to estimate the benefit of employing an eddy-resolving grid within the assimilation step, as compared with an eddy-permitting one, while keeping the eddy-resolving grid within the forecast step. It is shown that the use of a coarser grid for data assimilation does not allow an optimal exploitation of the present remote sensing observation network. A global decrease of about 15% in the error statistics is found when assimilating dense surface observations, and no significant improvement is seen for sparser observations (in situ profilers).