I'v started to look closer at curve25519, and I've added support for it in the eccdata program.
For the ecc operations, my current thinking is that one should use the Edwards curve equivalence outlined in http://cr.yp.to/papers.html#newelliptic, rather than the Montgomery x-coordinate method suggested in the curve25519 paper. The x-coordinate method needs an addition chain with additional values, which is a bit alien to all other scalar ecc multiplication in Nettle. The equivalent Edwards curve is
x^2 + y^2 = 1 + (121665/121666) x^2 y^2 (mod p).
Computations should be about as fast (according to the paper), and since the constant (121665/121666) is a non-square, the formulas are "complete", with no need to handle any special cases.
One needs conversion of the coordinates, and one also needs a square root to get the y coordinate, but I hope that shouldn't be too difficult.
Regards, /Niels