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.