The two-dimensional harmonic polylogarithms G(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=<(y+z)=<1. This algorithm is implemented into a FORTRAN subroutine tdhpl to compute two-dimensional harmonic polylogarithms up to weight<space>4.