We perform a nonperturbative determination of the O(a)-improvement coefficient c(SW) 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 8(3) X 16 employing the Schrodinger functional setup of lattice QCD. As Our main result, we obtain an interpolation formula for c(SW) and the critical hopping parameter K-c 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 c(SW) suggests that finite volume effects in c(SW) are not large on the 8(3) X 16 lattice. We investigate N-f dependence of c(SW) by additional simulations for N-f = 4, 2, and 0 at beta = 9.6. As a preparatory step for this study, we also determine c(SW) 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 c(SW) has been performed.