Support NEMO tripolar T and F pivots. Add option to follow latitude circles

[geoffstanley]

New options for the topology argument:

• topology = "tripolar-Tpivot" handles the NEMO tripolar grid pivoting on the T (tracer) point
• topology = "tripolar-Fpivot" handles the NEMO tripolar grid pivoting on the F (vorticity) point

The formulas for how to cross the north fold were created by studying Appendix E of the NEMO manual, v4.2.

New argument curve:

• curve = "great circle" is the default; no changes from previous behaviour
• curve = "latitude circle" is new; the walk through grid space will approximate a latitude circle.

Minor updates to examples/4_sections_on_global_tripolar_grid.ipynb to showcase curve = "latitude circle".

Minor edits to docs/environment.yml and examples/load_example_model_grid.py that I needed to make to get started.

[geoffstanley]

I left topology = "MOM-tripolar" unchanged for backwards compatibility. Is it correct? Should it be the same as one of tripolar-Tpivot or tripolar-Fpivot?

I found the MOM5 elements manual section 9.4.3 specifies that the northern seam for U points (same as F points, as that document is written for an Arakawa B-grid) operates as
U_{i, j} = U_{ni - i, 2 nj - j} for nj < j <= nj + halo
in 1-indexing. Plugging in j = nj + 1 gives
U_{i, nj+1} = U_{ni - i, nj - 1}
and then converting to 0-indexing, I think this becomes
U_{i, nj} = U_{(ni - 2 -i) % ni, nj-2}
which is the same formula I have obtained from the NEMO manual linked above. Is MOM6 the same as MOM5 in this regard?

From examples/4_sections_on_global_tripolar_grid.ipynb, here is the sections that was noted as causing trouble for Sectionate, with "MOM-tripolar":

topology	results
"MOM-tripolar"
"tripolar-Tpivot"
"tripolar-Fpivot"

Oddly, tripolar-Tpivot appears to work well, but MOM-tripolar and tripolar-Fpivot do not.

I thought the issue might be with the (i,j) walk algorithm, so I made a separate branch with a new algorithm over in https://github.com/geoffstanley/sectionate/tree/dev--walk-algorithm . More on that later.

I get good results with following various transects across the north fold with a mesh_mask from NEMO's eorca1 ("tripolar-Fpivot") and eorca025 ("tripolar-Tpivot") grids.

[hdrake]

Thanks for these great contributions, @geoffstanley !

Reviewing your PR inspired me to finally put in the effort (assisted by @claude) to do a major overhaul of sectionate's internals so it finds its neighbors using the grid-topology-aware internals of xgcm as opposed to reinventing the wheel for this inside of sectionate. While your patch here seems to work for several tripolar grids, it is hard to see how we would ever support more complicated grids like a cubed-sphere or LLC.

The major overhaul PR for sectionate is #47, which only extends support to the tripolar grid thanks to the new xgcm PR xgcm/xgcm#711, which works by extending the internal grid._pad method so that it correctly computes the halo padding across the tripolar grid seam (for both scalars and vectors).

I am going to hold off on this PR until I get further test and get feedback on the other two PRs, since if they are merged they would supercede this one.

[geoffstanley]

Ok, that sounds good. It certainly makes sense to have all the grid connections handled by xgcm. I see you implemented curve="latitude circle" so that's great too.

Regarding my other branch dev--walk-algorithm, as I recall from my testing of CMIP6/7 sections (eg for mfo and simassacrossline) in eorca1 and eorca025 NEMO grids, the new1 algorithm gave identical results to your original algorithm. Of course that doesn't mean it will give identical results in all cases, so more testing would be required. If they actually are identical, it would just be to improve the clarity of the algorithm (in my opinion, but I think you would agree!). If they aren't, then with a case in hand we could see which one is "correct" or better. I'll leave it to you to decide if you're interested in it; plus I think you have a larger collection of tests, the better to detect any difference between the algorithms. If you decide you want a PR for it, just let me know.

[hdrake]

@geoffstanley , I've done some more testing of my own implementation of the latitude-circle pathing. Can you explain what your intended behavior is if the user uses the curve='latitude-circle' option for a meridional segment or even just any section that is not purely zonal? Should the latitude-circle option only walk along constant latitude lines if the target points are at the same latitude, and otherwise default to the great-circle behavior?

[hdrake]

Related comment: https://github.com/MOM6-community/sectionate/pull/47/changes#r3752736092

[geoffstanley]

Hey @hdrake, my intention was to throw an error if the user requests curve='latitude-circle' and the endpoints of the segment have different latitudes (a purely zonal segment). But sorry, it looks like I didn't actually code that raise in my PR here.

Should the latitude-circle option only walk along constant latitude lines if the target points are at the same latitude, and otherwise default to the great-circle behavior?

Yes, I'd say so.

It does not make sense to use the latitude-circle to do a meridional segment. Just use a great circle between the same endpoints for that.

[hdrake]

@geoffstanley , I am thinking of testing this PR and potentially releasing it as 0.3.5, which leaves open the option that a future 0.4.0 release would extend support to multi-tile grids like ECCOv4r4 (e.g. as in my PR#47 that leverages new xgcm features for finding neighbors across face-connections and tripolar folds).

Before merging this, I'd like to add the following two changes:

1. Raise error if a segment has two endpoints with different latitudes under the "latitude-circle" option
2. Add the "latitude-and-great-circle" option I mentioned in the other PR

What do you think of this plan?

[geoffstanley]

That 0.3.5 plan sounds good to me.

1. I can add this (raise an error) to this PR.
2. I can probably do that too, just let me know if you do want it based on the other discussion.

[geoffstanley]

I've done 1. Let me know about 2.! And any other changes that would help make executing the PR more seamless. Eg. "great circle" or "great-circle"? I did the former but a lot of the discussion has used the latter.

[hdrake]

I think stick we should stick with the current API of space-separated e.g. "great circle". I had just misremembered.

I think it is worthwhile to add (2). Your suggested workaround of using join_sections also works but I think "latitude and great circle" (open to better names!) will be popular enough for users to add it as a supported option. I also can imagine users wanting to make sections with a long list of target coordinates that would be annoying to manually split up and recombine using join_sections.

[hdrake]

ReviewNB wasn't working but here are two comments about the notebook changes:

[geoffstanley]

I had read up on the MOM5 tripolar grid (I could not find a MOM6 description, so I assumed it was the same). From the MOM5 manual, I gleaned that the up neighbour of an F point at the bifold worked like this:
up = (2*nj-j, ni-i) # 1-index based
Translating that to be 0-index based, I got
up = (j-1, (nx-2-i)%nx)
which is the same as what I had already done for the (NEMO) tripolar F-point pivot.
So I assumed that the MOM6 tripolar grid was an F-point pivot, the same as NEMO has.
See these comments in my other branch: dev-walk-algorithm section.py#L598

But, as I recall, using 'tripolar-Tpivot' seemed to work better in that case. I was experimenting, and it led me to creating the dev-walk-algorithm branch to try changing the fundamental grid walking algorithm. But in the end, for the sections I was testing, my new algorithm gave identical results as your original algorithm. Though I did think my alternate algorithm made a bit more sense, and I wonder if it could lead to different results in other cases.

I'll review what's in the notebooks on this branch.

[geoffstanley]

In examples/4_sections_on_global_tripolar_grid.ipynb, I've removed the tests of tripolar-Tpivot and added an example of curve='"latitude and great circle". The latter prompted me to move the logic for that to infer_grid_path_from_geo so that it tests the input latitudes, not the latitudes after moving them to the nearest point on the grid.

[geoffstanley]

Just a minor update: I now average the latitudes of the two grid points that define the ends of the curve to be followed in infer_grid_path. We only allow curve="latitude circle" when the specified latitudes at either end of the section are exactly equal, but in infer_grid_path the geographic endpoints are converted to their nearest grid points, which may not have exactly equal latitudes. Rather than comparing against just one of those latitudes, I now compare against their average.

I thought about instead passing the exact latitude of the section (latstart) to infer_grid_path, but decided against this since it currently doesn't pass the exact lat & lon to infer_grid_path when using "great circle" either.

(force push was just to fix execution order diffs in the notebook, sorry.)

I think it's ready to go now, if you want. I also ran pytest on it, all 13 tests passed

=================================== 13 passed, 10 warnings in 10.20s ===================================