  Inclusion of a term $-\gamma\nabla\nabla\cdot u$, forcing $\nabla\cdot u$ to
be pointwise small, is an effective tool for improving mass conservation in
discretizations of incompressible flows. However, the added grad-div term
couples all velocity components, decreases sparsity and increases the condition
number in the linear systems that must be solved every time step. To address
these three issues various sparse grad-div regularizations and a modular
grad-div method have been developed. We develop and analyze herein a synthesis
of a fully decoupled, parallel sparse grad-div method of Guermond and Minev
with the modular grad-div method. Let
$G^{\ast}=-diag(\partial_{x}^{2},\partial_{y}^{2},\partial_{z}^{2})$ denote the
diagonal of $G=-\nabla\nabla\cdot$, and $\alpha\geq0$ an adjustable parameter.
The 2-step method considered is $$\begin{eqnarray} 1
&:&\frac{\widetilde{u}^{n+1}-u^{n}}{k}+u^{n}\cdot \nabla
\widetilde{u}^{n+1}+\nabla p^{n+1}-\nu \Delta \widetilde{u}^{n+1}=f\text{ &
}\nabla \cdot \widetilde{u}^{n+1}=0,\\ 2 &:&\left[ \frac{1}{k}I+(\gamma +\alpha
)G^{\ast }\right] u^{n+1}=\frac{1}{k }\widetilde{u}^{n+1}+\left[ (\gamma
+\alpha )G^{\ast }-\gamma G\right] u^{n}. \end{eqnarray}$$ We prove its
unconditional, nonlinear, long time stability in $3d$ for
$\alpha\geq0.5\gamma$. The analysis also establishes that the method controls
the persistent size of $||\nabla\cdot u||$ in general and controls the
transients in $||\nabla\cdot u||$ when $u(x,0)=0$ and $f(x,t)\neq0$ provided
$\alpha>0.5\gamma$. Consistent numerical tests are presented.
