  A rigorous generation of spin-adapted (spin-free) substitution operators for
high spin ($S=S_z$) references of arbitrary substitution order and spin quantum
number $S$ is presented. The generated operators lead to linearly independent
but non-orthogonal CSFs when applied to the reference and span the complete
spin space. To incorporate spin completeness, spectating substitutions (as e.g.
$\hat{E}_{iv}^{va}$) are introduced. The presented procedure utilizes
L\"owdin's projection operator method of spin eigenfunction generation to
ensure spin completeness. The generated operators are explicitly checked for
(i) their linear independence and (ii) their spin completeness for up to
10-fold substitutions and up to a multiplicity of $2S+1 = 11$. A proof of
concept implementation utilizing the generated operators in a coupled cluster
(CC) calculation was successfully applied to the high spin states of the Boron
atom. The results show pure spin states as well as small effects on the
correlation energy compared to spinorbital CC. A comparison to spin-adapted but
spin-incomplete CC shows a significant spin incompleteness error.
