When using all the information optimally, it's probably desirable to model the 1-D sky spectra for an arm (generating a 2-D model and making use of multiple dithered realisations to recover the sampling), subtract the sky, and then merge the 1-D sky spectra across the arms. The sky subtraction should be done pre-merge as the science spectra are not very well samples, and sky subtraction removes a lot of high-frequency power.
However, while we are using a simpler algorithm, merging the arms to create a per-arm sky then interpolating it to the per-fibre wavelength solutions means that we are interpolating the sky twice, and so with e.g. linear interpolation the sky lines have an extra sqrt(1/6) added in quadrature to the line widths. The simplest way to avoid this is to merge all the spectra (including the sky spectra) and then merge the sky spectra, and then subtract.

Please add an option doSubtractSky1dBeforeMerge to control this; for now the default should be False.