  We present an extension of the RSVD-$\Delta t$ algorithm initially developed
for resolvent analysis of statistically stationary flows to handle harmonic
resolvent analysis of time-periodic flows. The harmonic resolvent operator, as
proposed by \citet{Padovanetal20}, characterizes the linearized dynamics of
time-periodic flows in the frequency domain, and its singular value
decomposition reveals forcing and response modes with optimal energetic gain.
However, computing harmonic resolvent modes poses challenges due to $(i)$ the
coupling of all $N_{\omega}$ retained frequencies into a single harmonic
resolvent operator and $(ii)$ the singularity or near-singularity of the
operator, making harmonic resolvent analysis considerably more computationally
expensive than a standard resolvent analysis. To overcome these challenges, the
RSVD-$\Delta t$ algorithm leverages time stepping of the underlying
time-periodic linearized Navier-Stokes operator, which is $N_{\omega}$ times
smaller than the harmonic resolvent operator, to compute the action of the
harmonic resolvent operator. We develop strategies to minimize the algorithm's
CPU and memory consumption, and our results demonstrate that these costs scale
linearly with the problem dimension. We validate the RSVD-$\Delta t$ algorithm
by computing modes for a periodically varying Ginzburg-Landau equation and
demonstrate its performance using the flow over an airfoil.
