  The two-dimensional harmonic polylogarithms $\G(\vec{a}(z);y)$, a
generalization of the harmonic polylogarithms, themselves a generalization of
Nielsen's polylogarithms, appear in analytic calculations of multi-loop
radiative corrections in quantum field theory. We present an algorithm for the
numerical evaluation of two-dimensional harmonic polylogarithms, with the two
arguments $y,z$ varying in the triangle $0\le y \le 1$, $ 0\le z \le 1$, $\
0\le (y+z) \le 1$. This algorithm is implemented into a {\tt FORTRAN}
subroutine {\tt tdhpl} to compute two-dimensional harmonic polylogarithms up to
weight 4.
