Mathematica/3DContours

Mathematica supports 3D contours on surfaces using the option MeshFunctions.

MeshFunctions draws ISO lines of an arbitary function of the x,y,z values, and for some functions additional parameters.

For a conventional set of altitude contours we just use the identity function of the z value (the third parameter) using MeshFunctions->{#3&}

eg

Plot3D[Sin[x y], {x, 0, 3}, {y, 0, 3}, MeshFunctions -> {#3 &}]

In addition the region between the contours can be colored using the option MeshShading. The following example divides the image into 10 altitude contours and colors them in gray scale.

Plot3D[Sin[x y], {x, 0, 3}, {y, 0, 3}, Mesh -> 10,  MeshFunctions -> {#3 &},  MeshShading -> GrayLevel /@ Range[0, 1, 0.1], Lighting -> "Neutral"]