phase wrapping group delay
The relationship between group delay and phase is

\$t_g(\omega) = -\dfrac{d\phi}{d\omega}\$.

So to get phase from group delay you don't subtract, you integrate.
Since these are all real numbers (not complex), you can do this integration numerically by any standard method, such as trapezoidal or Simpson's rule. If you have access to a numerical math package like Mathematica, Matlab, or Octave, it will have a built in routine to do this integration.

Once you do this integration numerically, there's nothing that says the phase you extract will be banded between 0 and 2\$\pi\$. So there's no reason you should have to unwrap the phase. 

If you do the measurement the other way around (that is, measure phase and then calculate phase delay and group delay), then you typically measure phase between 0 and 2\$\pi\$ and you have to unwrap it to obtain a group delay without glitches where the raw phase measurement is discontinuous.