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≤y≤1, 0≤z≤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.Comment: 22 pages, LaTe