We non-perturbatively determine the renormalization factor of the axial vector current in lattice QCD with $N_f=3$ flavors of Wilson-clover fermions and the tree-level Symanzik-improved gauge action. The (by now standard) renormalization condition is derived from the massive axial Ward identity and it is imposed among Schr{o}dinger functional states with large overlap on the lowest lying hadronic state in the pseudoscalar channel, in order to reduce kinematically enhanced cutoff effects. We explore a range of couplings relevant for simulations at lattice spacings of $approx 0.09$ fm and below. An interpolation formula for $Z_A(g_0^2)$, smoothly connecting the non-perturbative values to the 1-loop expression, is provided together with our final results.