We perform a nonperturbative determination of the O(a)-improvement coefficient cSW for the Wilson quark action in three-flavor QCD with the plaquette gauge action. Numerical simulations are carried out in a range of beta=12.0–5.2 on a single lattice size of 83×16 employing the Schrödinger functional setup of lattice QCD. As our main result, we obtain an interpolation formula for cSW and the critical hopping parameter Kc as a function of the bare coupling. This enables us to remove the O(a) scaling violation from physical observables in future numerical simulation in the wide range of beta. Our analysis with a perturbatively modified improvement condition for cSW suggests that finite volume effects in cSW are not large on the 83×16 lattice. We investigate Nf dependence of cSW by additional simulations for Nf=4, 2, and 0 at beta=9.6. As a preparatory step for this study, we also determine cSW in two-flavor QCD at beta=5.2. At this beta, several groups have carried out large-scale calculations of the hadron spectrum, while no systematic determination of cSW has been performed.