cosima-recipes
cosima-recipes copied to clipboard
xgcm vorticity formula
@navidcy should https://github.com/COSIMA/cosima-recipes/blob/master/DocumentedExamples/RelativeVorticity.ipynb
use dxt
instead of dxu
in cell 18?
ζ_xgcm = ( grid.interp( grid.diff(ds.v, 'X') / grid.interp(ds.dxu, 'X'), 'Y', boundary='extend')
- grid.interp( grid.diff(ds.u, 'Y', boundary='extend') / grid.interp(ds.dyt, 'X'), 'X') )
dxu
is inconsistent with the subsequent text.
What do you mean "inconsistent with the subsequent text"?
The objective in this cell is to reproduce ζ_mom5
. If you can do that with dyu
then go for it and make a PR.
I'm not suggesting that what I did is definitely the only way forward. I don't even recall my thinking. But at least by plotting the difference you can see that (the values at least) of ζ_xgcm
match those of ζ_mom5
. The plot shows almost-zero (see the figure below). However, ζ_xgcm
lives on t-grid while ζ_mom5
lives on u-grid.
I haven't looked at this in great detail. It just seems odd to have dxu
and dyt
in the same formula, especially as the "simpler formula" given afterwards uses dxt
and dyt
:
ζ_xgcm = grid.interp(grid.diff(ds.v, 'X'), 'Y', boundary='extend')/ds.dxt - grid.interp(grid.diff(ds.u, 'Y', boundary='extend'), 'X')/ds.dyt
It's odd. If you can help sort it out in a cleaner way it'll be great.
But the objective of the first (seemingly weird) method was to replicate was MOM5 is actually doing. I couldn't figure out a way to do that using dyt
! And also I couldn't figure how to do this computation and end up on the u-grid (like where ζ_mom5
lives). That is, the way the notebook is now ζ_mom5
and ζ_xgcm
have the same values but live on different grids! If you can figure out how to replicate the values of ζ_mom5
and end up on the u-grid let me know. However, what (I think I understood that) you suggested doesn't even give same values!
ζ_xgcm.values[200:205, 100:104]
array([[ 1.1507019e-06, 1.0019297e-06, 5.4093033e-07, 2.8089801e-07],
[ 9.9383203e-07, 9.8258420e-07, 1.0025971e-06, 9.3004383e-07],
[ 1.9963846e-07, 5.3464299e-07, 7.5498042e-07, 4.4450223e-07],
[-1.6323114e-07, -1.6445910e-08, 2.2685566e-07, -1.9235793e-07],
[-6.1921597e-08, 2.8606735e-08, 3.5336069e-07, 2.6428950e-08]],
dtype=float32)
ζ_aekiss = (grid.interp( grid.diff(ds.v, 'X') / grid.interp(ds.dxu, 'X'), 'X')
- grid.interp( grid.diff(ds.u, 'Y', boundary='extend') / grid.interp(ds.dyu, 'Y', boundary='extend'), 'Y', boundary='extend'))
-
ζ_aekiss.values[200:205, 100:104]
array([[ 1.0400078e-06, 8.1025718e-07, 6.4869403e-07, 9.7753821e-07],
[ 7.0336637e-07, 8.3177400e-07, 8.4166714e-07, 7.4260828e-07],
[ 1.4705354e-07, 4.5017367e-07, 3.7383137e-07, -1.6373960e-07],
[-4.7659199e-08, 1.5902492e-07, 1.2845433e-07, -4.5021255e-07],
[-3.4233143e-08, 2.7488650e-07, 3.5775838e-07, -1.7128500e-08]],
dtype=float32)
ζ_mom5.values[200:205, 100:104]
array([[ 1.1507033e-06, 1.0019312e-06, 5.4093209e-07, 2.8090085e-07],
[ 9.9383305e-07, 9.8258579e-07, 1.0025997e-06, 9.3004758e-07],
[ 1.9963758e-07, 5.3464328e-07, 7.5498144e-07, 4.4450331e-07],
[-1.6323250e-07, -1.6447160e-08, 2.2685441e-07, -1.9235949e-07],
[-6.1921909e-08, 2.8606166e-08, 3.5335984e-07, 2.6427927e-08]],
dtype=float32)
But if replicating the MOM5 diagnostic is not part of the requirements then you should just use the "simpler"
ζ_simpler = grid.interp(grid.diff(ds.v, 'X'), 'Y', boundary='extend')/ds.dxt - grid.interp(grid.diff(ds.u, 'Y', boundary='extend'), 'X')/ds.dyt
which lives on the t-grid, or the even simpler
ζ_evensimpler = ds.v / ds.dxu - ds.u / ds.dyu
which lives on the u-grid.
Bottom line is that it all depends what you wanna compute. Maybe I have missed your point entirely? Sorry if I have done that... But changing dxt
->dxu
just for the sake of symmetry ignoring
I think you have a point here. Or at least there is something fishy.... because the ζ_xgcm
lives on the cell centers while ζ_mom5
lives at the corners.
Question: the vorticity_z
diagnostic from MOM5 in which location is it output??
Is this still relevant @aekiss ?