Re: [ferret_users] potential vorticity from OGCM

[Ryo Furue <furue@xxxxxxxxxx>]

Hi Prescilla,

| The code I prepared is as below
| 
| use 2014_hori_U_AM500_AH250.nc
| use 2014_hori_V_AM500_AH250.nc
| set region/X=78E:94E/Y=6N:23N/L=8/K=10
| use density.nc
| LET zeta=(VVEL[d=2,x=@ddc])-(UVEL[d=1,y=@ddc])  ! relative vorticity
| let deg2rad=(2*3.14)/360
| let omega = 7.292e-5
| let f=2*omega*sin(deg2rad*Y[g=UVEL[d=1]]) !planetary V
| let rho_zero=1026 ! referece density
| let grad=DENSITY[z=@ddc] ! density gradient
| fill (grad*(zeta+f))/rho_zero

Looks good to me, if DENSITY, VVEL, and UVEL are in MKS units.
A potential (minor) problem is that if DENSITY, VVEL, and UVEL
are defined on different grids, then it's better to map them on
to a common grid before the final multiplication; for example

  let pv = grad[g=zeta] * (zeta + f) / rho_zero

Density should be a potential density but I don't think
that matters much for the upper ocean.

By the way, I have a very minor piece of advice, which isn't
related to the cause of your problem.  I think it's better
not to throw away precision gratuitously: I would write
something along the lines of

  let pi = 4*atan(1.0)  !! or acos(-1.0)
  let deg2rad = pi/180
  let omega = 2*pi/(24*60*60)

Regards,
Ryo