The Multi-region Relaxed MHD (MRxMHD) has been successful in the construction of equilibria in three-dimensional (3D) configurations. In MRxMHD, the plasma is sliced into sub-volumes separated by ideal interfaces, each undergoing relaxation, allowing the formation of islands and chaos. The resulting equilibrium has a stepped pressure profile across sub-volumes. The Stepped Pressure Equilibrium Code (SPEC) [S.R. Hudson et al , Phys. Plasmas 19, 112502 (2012)] was developed to calculate MRxMHD equilibria numerically. In this work, we have extended the SPEC code to compute MRxMHD equilibria with field-aligned flow and rotation, following the theoretical development to incorporate cross-helicity and angular momentum constraints. The code has been verified for convergence and compared to a Grad-Shafranov solver in 2D. We apply our new tool to study the flow profile change before and after the sawtooth crash of a reversed-field pinch discharge, in which data of the parallel flow is available. We find the promising result that under the constraints of cross-helicity and angular momentum, the parallel flow profile in postcrash SPEC equilibrium is flat in the plasma core and the amplitude of the flow matches experimental observations. Finally, we provide an example equilibrium with a 3D helical field structure as the favoured lower energy state. This will be the first 3D numerical equilibrium in which the flow effects are self-consistently calculated.