Hi Rachel,
Sorry for the lack of response – it’s always hard to find the headspace to think and provide adequate answers. Hopefully this will be useful anyway.
Since that original post 5 years ago, @rsmith has introduced the dwicat command / script, which is expressly designed for this use case, taking into account any differences in scaling that might be present between the acquisitions due to e.g. scanner recalibration. It’ll accept all four of your scans in the same command if that’s what you want to do – more on that below.
Regarding the subsequent stages in the pipeline, the Gibbs ringing removal operates on a slice-wise basis and won’t be affected at all however your data are combined.
However, the denoising operates differently, given that is looks for redundancy in the entire signal series over neighbouring voxels. There have been a few threads looking at whether it’s best (or even just OK) to concatenate different PE directions prior to denoising, or whether it’s best to denoise each series separately (e.g. this thread), though I can’t find any that address the case where the full DW series is acquired for each PE direction (if anyone can find such posts, please link to them here!). Personally, I reckon I would concatenate the lot and denoise that, but you may want to double-check what difference that really makes (e.g. by comparing the RMS differences between the full series pre & post denoising, in one case with each PE direction denoised separately then combined, in the other with all PE directions combined prior to denoising). If you go to that trouble, I (and no doubt many others) would be very interested to hear about your findings… 
The last stage would naturally be dwifslpreproc, and here the question of the order in which you concatenate your data becomes important. Your acquisition fits into the -rep_all scenario (described in this section of the documentation for dwifslpreproc – look under the heading “All DWI directions & b-values are acquired twice, with the phase encoding direction of the second acquisition protocol being reversed with respect to the first”). This option assumes that your input DWI data consists of two repeats, with the first half acquired with one PE direction, the second half with the opposite PE direction. As long as your data have been concatenated accordingly, this should work out of the box.
In summary, I reckon you could process your data using a pipeline like this:
dwicat subjid_91directions_RL.mif subjid_90directions_RL.mif subjid_91directions_LR.mif subjid_90directions_LR.mif all_dwi.mif
dwidenoise all_dwi.mif all_dwi_dn.mif -noise noise_map.mif
mrdegibbs all_dwi_dn.mif all_dwi_dn_dg.mif
dwifslpreproc all_dwi_dn_dg.mif dwi_out.mif -rpe_all -pe_dir rl -readout_time 0.66
You may wish to edit the readout time to match your acquisition, though I don’t think this will make any difference to the outcome in your case.
This is obviously based on my gut feeling that it’s OK to denoise the full concatenated set, so feel free to experiment and amend as required!
All the best,
Donald.