Intel Fortran has a compiler switch
to treat all double precision
variables as 128 bit. It is

-double_size 128

All I did was to replace the
"real*8" statements by
"double precision" in the *.f and *.h
and recompile with Intel Fortran.
It is necessary to reduce the
predefined machine epsilon for the
NAG library routines. It is in
x02ajft.f. I simply changed one line
in this file to

X02AJF = 1.e-22

Normally, one would expect that a mantissa
of more than 30 decimal places is supported,
but this turned out not to be the case
on a Pentium M and Intel Fortran V.8
(I experimented a bit to determine
the effective length of the mantissa).
However, on an Athlon64 with Intel Fortran 9.1
the mantissa seems to have the 30+ decimal
places one would expect for 128bit.