I was under the impression that if you took some 3D co-ordinates, and fed
them into screenX/Y/Z that the resulting numbers would give you the
original 3D co-ordinates if fed back into modelX/Y/Z however this does not
seem to be the case.

e.g.:
void setup()
{
size(500,400,P3D);

float sX=screenX(200,200,-200);
float sY=screenY(200,200,-200);
float sZ=screenZ(200,200,-200);
float mX=modelX(sX,sY,sZ);
float mY=modelY(sX,sY,sZ);
float mZ=modelZ(sX,sY,sZ);
println("Screen ("+sX+","+sY+","+sZ+") Model ("+mX+","+mY+","+mZ+")");
}

gives an output of:
Screen (218.30127,200.0,0.94606316) Model (-31.69873,0.0,-345.46408)

However the Model co-ordinates shoud be (200,200,-200)

The 0.0 for the Y value looks interesting, as if when converting
screen->model a translate(width/2,height/2); has been assumed, which would
also partly explain wht the X is -31.69. However that doesnt' explain the
difference in the modelZ. Possibyl the processing internal z-translate (so
that widht/2,height/2,0 is not right on top of the camera, but some
distance away)

no, you're misunderstanding their intent. they just give you the transformed versions of the
points at various stages in the rendering pipeline.

Be as it may, we certainly would like to have the dual functions, sort of
what gluProject and gluUnproject do.