Better chunking strategies for constlat intersections and zonal routines.

[cmdupuis3]

This PR contains two post-accusphere optimizations, eliminating low-level hard materializations by using a vector-based masking strategy rather than individual conditionals, and reduced zonal_mean peakmem by building only the candidate faces instead of the whole grid.

Partly addresses #1587

Waiting on accusphere merge to main #1513 first.

Overview

Phase A eliminated the largest peak-memory amplifier in the zonal path: the port materialized the whole-grid (n_face, n_max, 2, 3) face-edge array up front (a 5.8× blow-up over the node coordinates, with a ~23 MB build peak on a 28k-face grid) even though each latitude only touches ~1% of faces. I added an @njit(parallel=True) subset builder, _get_cartesian_face_edge_nodes_array_subset, that constructs edges only for the candidate faces of a given latitude/band — bit-identical to indexing the full array — and rewired _compute_non_conservative_zonal_mean and _compute_face_band_weights to build per-candidate subsets instead of the whole grid. Peak memory for a 180-latitude zonal_mean dropped from 23.66 MB to 0.27 MB (≈88×) and it ran ~10% faster (no full build, no per-latitude fancy-index copies), verified lossless via a git-stash A/B (the only diff was a pre-existing 4.4e-16 parallel-reduction nondeterminism) with the full zonal/cross-section suite passing.

Phase B moved the six edge/face screeners (constant_lat/lon_intersections_no_extreme, constant_lat/lon_intersections_face_bounds, faces_within_lat/lon_bounds) off @njit and onto plain vectorized NumPy, drawing the boundary so the low-level Numba kernels stay for real geometry while these memory-bound elementwise predicates use NumPy — which is ~2.1× faster here and, unlike an njit kernel that forces a full .values materialization, composes with dask (a single _flatnonzero helper leans on NumPy's array_function protocol so a dask mask reduces block-wise, no explicit dask branch). Call sites pass .data instead of .values so a chunked grid stays lazy, and edge coordinates are gathered positionally (node_z.data[conn.ravel()].reshape(...)) to stay chunk-friendly. Along the way it fixed two real latent bugs: get_edges_at_constant_latitude referenced a nonexistent self.edge_node_z (it raised AttributeError on every call), and both edge paths crashed on chunked grids because xarray can't vindex with a dask indexer. Results are bit-identical to the original per-element loops across 300 randomized trials, dask==numpy, and it's committed as bf0abbe "Lazy intersections".

PR Checklist

General

• An issue is linked created and linked
• Add appropriate labels
• Filled out Overview and Expected Usage (if applicable) sections

Testing

• Adequate tests are created if there is new functionality
• Tests cover all possible logical paths in your function
• Tests are not too basic (such as simply calling a function and nothing else)

Documentation

• Docstrings have been added to all new functions
• Docstrings have updated with any function changes