A two-step approximate factorization technique for implementing a computationally stable nonlinear unsteady frequency-domain harmonic balance solution method within existing implicit computational fluid dynamic flow solver codes is presented. The approach uses an explicit discretization of the harmonic balance source term, and no new implicit code development is required. Both of these features enable the harmonic balance method to be implemented within existing implicit flow solver codes with minimal modification necessary to the underlying flow solver code. The resulting harmonic balance solver can then be used for modeling nonlinear periodic unsteady flows. The methodology is applied to the NASA OVERFLOW flow solver code, and results are presented for transonic viscous flow past an unsteady pitching airfoil section.