Interpolating a scalar field in a 3D space
Since you have a spatial mesh with constant spacing, you can identify all neighbors on opposite sides of the isosurface.  Choose some form of interpolation (q.v. Reed Copsey's answer) and do root-finding along the line between each such neighbor.
I have a 3D space (x, y, z) with an additional parameter at each point (energy), giving 4 dimensions of data in total.
I would like to find a set of x, y, z points which correspond to an iso-energy surface found by interpolating between the known points.
The spacial mesh has constant spacing and surrounds the iso-energy surface entirely, however, it does not occupy a cubic space (the mesh occupies a roughly cylindrical space)
Speed is not crucial, I can leave this number crunching for a while. Although I'm coding in Python and NumPy, I can write portions of the code in FORTRAN. I can also wrap existing C/C++/FORTRAN libraries for use in the scripts, if such libraries exist.
All examples and algorithms that I have so far found online (and in Numerical Recipes) stop short of 4D data.