Calculate geostrophic velocities from sea level height values

[florenciab]

I have a climatological year's average of ocean variables obtained from a numerical model. I would like to use the sea level height values to calculate geostrophic velocities (ug, vg) from the balance between Coriolis Force and Pressure Gradient Force in RStudio. I have already calculated the Coriolis parameter (f) and I have de g value, but I am stuck in calculating the sea level height variations in x (dz/dx) and y (dz/dy). I understand that the "grad" function of the "numDeriv" package could help me, but I don't know how to use it. Any suggestions?

x <- lon * 60 * 1852 * cos(lat*pi/180)  
y <- lat * 60 * 1852
omega <- 2*pi/(24*60*60)    
g <- -9.8
f <- 2 * omega * sin(lat*pi/180)
ug <- -(g/f)*(dz/dy)
vg <- (g/f)*(dz/dx)

I am transcribing a Matlab code, where they calculate the velocities using the function "nderiv_fornberg" , but I have not found equivalents in R. I leave below the matlab code.

[si, sj]=size(ssh);
for j=1:sj
dzeta_dx(:,j)=nderiv_fornberg( 1, x(:,j), ssh(:,j));
end
for i=1:si
dzeta_dy(i,:)=nderiv_fornberg( 1, y(i,:), ssh(i,:));
end

where "ssh" is the matrix containing the values of z (sea surface height).

Thanks in advance.