Domain decomposition procedures for solving parabolic equations are considered. The underlying discretization is block-centered finite differences. In this procedures, fluxes at subdomain interfaces are calculated from the solution at the previous time level. These fluxes serve as Neumann boundary data for implicit, block-centered discretizations in the subdomains. A priori error estimates are derived, and numerical results examining the stability, accuracy, and parallelism of the scheme are presented.