The coma of comet C/2016 R2 (PanSTARRS) is one of the most chemically peculiar ever observed, in particular due to its extremely high CO/H2O and N2+ /H2O ratios, and unusual trace volatile abundances. However, the complex shape of its CO emission lines, as well as uncertainties in the coma structure and excitation, has lead to ambiguities in the total CO production rate. We performed high-resolution, spatially, spectrally, and temporally resolved CO observations using the James Clerk Maxwell Telescope and Submillimeter Array to elucidate the outgassing behavior of C/2016 R2. Results are analyzed using a new, time-dependent, three-dimensional radiative transfer code (SUBlimating gases in LIME; SUBLIME, based on the open-source version of the LIne Modeling Engine), incorporating for the first time, accurate state-to-state collisional rate coefficients for the CO-CO system. The total CO production rate was found to be in the range of (3.8 - 7.6) × 1028 s-1 between 2018 January 13 and February 1 (at r H = 2.8-2.9 au), with a mean value of (5.3 ± 0.6) × 1028 s-1. The emission is concentrated in a near-sunward jet, with a half-opening angle of ∼62° and an outflow velocity of 0.51 ± 0.01 km s-1, compared to 0.25 ± 0.01 km s-1 in the ambient (and nightside) coma. Evidence was also found for an extended source of CO emission, possibly due to icy grain sublimation around 1.2 × 105 km from the nucleus. Based on the coma molecular abundances, we propose that the nucleus ices of C/2016 R2 can be divided into a rapidly sublimating apolar phase, rich in CO, CO2, N2, and CH3OH, and a predominantly frozen (or less abundant), polar phase containing more H2O, CH4, H2CO, and HCN.