In this article, a parallel Conjugate Gradient Squared (CGS) block-centered finite difference scheme is introduced and analyzed to cast about the numerical solution of a nonlinear time-fractional parabolic equation with the Neumann condition and a nonlinear reaction term. The unconditionally stable result, which just depends on initial value and source item, is derived. Some a priori estimates of discrete L-2-norm with optimal order of convergence O(Delta t(2-alpha) + h(2) + k(2) ) with pressure and velocity are established on both uniform and non-uniform rectangular grids, where Delta t is the time step, h and k are maximal mesh sizes of x and y-directional grids. In our model, because the simulation is iterated over a series of time steps, it is most beneficial if all the calculation units in a time step can be run separately. In CGS algorithm, adding OpenMP instructions to the circulation calculations in an iterative step can realize parallel computing. To examine the efficiency and accuracy of the proposed method, numerical experiments using the schemes are studied. The results clearly show the benefit of using the proposed approach in terms of execution time reduction and speedup with respect to the sequential running in a single thread. (C) 2016 Elsevier B.V. All rights reserved.