This adapts the ticks and the mesh to the range of the $z$ values:
data2 = Table[10 + (3 ((j - 20)^2 - (j + 5)^2))/(75 + (i - 10)^2 + 2 (j - 10)^2),
{i, 0, 20}, {j, 0, 20}];
(* computes ticks and stores the coordinates in zmesh *)
ticks[zmin_, zmax_] := Module[{mark, ldz, dz},
ldz = Round[Log10[(zmax - zmin)/10], 1/3];
dz = Round[10^ldz, 10^Floor[ldz]];
mark = Round[10^(ldz + 2/3), 10^Floor[ldz + 2/3]];
If[Chop@Mod[#, mark] == 0, #, {#, "", {0.005, 0}}] & /@ (zmesh =
N@Range[Ceiling[zmin, dz], zmax, dz])];
myTicks = ticks[Min[data2], Max[data2]]; (* sets mesh & ticks to same z-values *)
ListPlot3D[data2, MeshFunctions -> {#3 &}, Mesh -> {zmesh},
BoxRatios -> {1, 1, 1/2}, PlotRangePadding -> 0, FaceGrids -> All,
Ticks -> {Automatic, Automatic, myTicks}]

This is kind of cute, except ListPlot3D
evaluates twice,
since zmesh
is reset by the function ticks
.
But perhaps that's ok if you have the spare time.
Dynamic@ListPlot3D[data2, MeshFunctions -> {#3 &}, Mesh -> {zmesh},
BoxRatios -> {1, 1, 1/2}, PlotRangePadding -> 0, FaceGrids -> All,
Ticks -> {Automatic, Automatic, ticks}]
It also works with other plot functions
Dynamic@Plot3D[5 Sin[x^2 + x y^2/2]/(1 + x^2 + y^2),
{x, -3, 3}, {y, -3, 3}, MaxRecursion -> 3, MeshFunctions -> {#3 &}, Mesh -> {zmesh},
BoxRatios -> {1, 1, 1/2}, PlotRangePadding -> 0, FaceGrids -> All,
Ticks -> {Automatic, Automatic, ticks}]
