vectorize 2D gradient with spatially varying bins
Your f and ssh are already vectorized, you can do the same quite trivially with u and v also. There is nothing tricky going on in your loop. The process is simply to remove the for, leaving the assignment i=2:size(x,1)-1. And replace all matrix multiplication and division by element-wise multiplication and division (.* and ./). This leaves:

%Calculate geostrophic current
u = zeros(size(ssh));
v = zeros(size(ssh)); 
i = 2:size(x,1)-1;
j = 2:size(y,2)-1;
dx(i,j) = (x(i+1,j)-x(i-1,j)) .* (R*cosd(y(i,j))*pi/180);
dy(i,j) = (y(i,j+1)-y(i,j-1)) .* (R*pi/180);
u(i,j) = -g(i,j) ./ f(i,j) .* (ssh(i,j+1)-ssh(i,j-1)) ./ dy(i,j);
v(i,j) =  g(i,j) ./ f(i,j) .* (ssh(i+1,j)-ssh(i-1,j)) ./ dx(i,j);


You can then do a slight simplification, dx and dy do not need indexing, since you're using the same part that you assign:

dx = (x(i+1,j)-x(i-1,j)) .* (R*cosd(y(i,j))*pi/180);
dy = (y(i,j+1)-y(i,j-1)) .* (R*pi/180);
u(i,j) = -g(i,j) ./ f(i,j) .* (ssh(i,j+1)-ssh(i,j-1)) ./ dy;
v(i,j) =  g(i,j) ./ f(i,j) .* (ssh(i+1,j)-ssh(i-1,j)) ./ dx;