@J.M I tried to implement your suggestion. I have the following problem. Starting from a 4x4x4 field, nn = 4;
u = GaussianRandomField[nn, 3, Function[k, Exp[-k^2/2]/k]] // Re;
uf = ListInterpolation[u];
fx = Function[{x, y, z}, D[uf[x, y, z], x] // Evaluate];
fy = Function[{x, y, z}, D[uf[x, y, z], y] // Evaluate];
fz = Function[{x, y, z}, D[uf[x, y, z], z] // Evaluate];
pl = ContourPlot3D[{fx[x, y, z], fy[x, y, z]}, {x, 1, nn}, {y, 1,
nn}, {z, 1, nn}, Contours -> {0}, ContourStyle -> Opacity[0],