I'm not sure how the coordinates in the gradient relate to longitude, but let's just say that the scaling parameters are xs,ys,zs and you define scalingParams = [xs ys zs]. Then you could create a function
function grad = gradScaled(f,scalingParams)
[gx,gy,gz] = gradient(f);
grad = [gx gy gz]./scalingParams;
g = @(f) gradScaled(f,s0);
and define
h = @(f) abs(g(f));
Then
tfp = dot(g(h(temp)),g(temp)/h(temp));
(note, by the way, that your first abs is in the wrong place).
