Retain XLONG_U/XLAT_V staggered edge in lon_m/lat_m

[jmsexton03]

This PR addresses https://github.com/erf-model/remora-erf-driver/issues/11

Specific changes / choices

• Grow typed subdomain so intersection copy keeps staggered edge
• Clamp lon_m/lat_m ghosts to var_fab extent, not cell domain
• Leave sinPhi_m/cosPhi_m on cell clamp; Coriolis reads their ghosts
• Document the staggering GetOceanToAtmosCornerCoordinates promises

Motivation Summary

XLONG_U and XLAT_V are edge-staggered in a wrfinput file — west_east_stag and south_north_stag, so nx+1 and ny+1 entries respectively. The fabs that receive them carry CELL index type, so the typed subdomain built at ERF_InitFromWRFInput.cpp:488 is a cell box, and the intersection copy that follows,

Box subdomain_to_fill_typed(convert(subdomain_tmp, var_fab_from_file.box().ixType()));
...
Box intersection = var_fab.box() & var_fab_from_file.box();

silently drops the final staggered column of XLONG_U and the final staggered row of XLAT_V. The ghost fill further down then clamps against the cell domain, so lon_m(nx,j) receives a duplicate of column nx-1 rather than the true eastern edge, and lat_m(i,ny) likewise.

Measured on a real wrfinput_d01 with west_east = 117, west_east_stag = 118:

	value
XLONG_U(0,50,117) in the file	-85.09668
lon_m(117,50) before this change	-85.38144 (duplicate of column 116)
error	0.285° — a full grid spacing

XLAT_V shows the same at the north edge (45.51647 in the file against 45.27907).

Why the outermost column matters

lon_m and lat_m are not only diagnostic. They are the coordinate source a coupled ocean model consumes as a cell-corner mesh, through ERF::GetOceanToAtmosCornerCoordinates. The consumer reconstructs true corners by averaging along the axis where each array is cell-centred and taking it as-is along the axis where it is already nodal:

// lon_m is already nodal in i and centered in j; lat_m the reverse.
node(i,j,0,0) = 0.5 * (lon(i,j-1,0) + lon(i,j,0));
node(i,j,0,1) = 0.5 * (lat(i-1,j,0) + lat(i,j,0));

That requires lon out to i = nx. With the duplicate in place, the easternmost column of corner quads collapses to zero width, its polygon area is zero, and every conservative-remap weight built against it vanishes — so the outermost column of the domain silently exchanges nothing. The north row fails the same way through lat_m.

The fix

Grow the typed subdomain by one in the staggered direction before the intersection copy, guarded so it can only ever widen into data the file actually contains:

if (var_name == "XLONG_U" &&
    var_fab_from_file.box().bigEnd(0) > subdomain_to_fill_typed.bigEnd(0)) {
    subdomain_to_fill_typed.growHi(0,1);
}

and clamp the lon_m / lat_m ghost fill against var_fab's own extent rather than the cell domain, so the retained column is what the ghost receives. var_fab keeps CELL index type; only its extent widens.

Note the apparently obvious fix — widening the clamp alone — is a no-op. Before this change var_fab.box().bigEnd(0) == i_hi == 116, so var_fab never contained column 117 and clamping against its own extent changed nothing. The subdomain has to grow first.

Why this approach, against established practice

This is worth stating because the alternative is what the reference implementation does, and the choice is not obvious.

SCRIP requires explicitly enumerated corners: a grid file carries grid_corner_lat(grid_size, grid_corners) and grid_corner_lon(...), four per cell for quads, boundary cells included (SCRIP 1.4 user guide; Jones, Mon. Wea. Rev. 127, 2204–2210, 1999, https://doi.org/10.1175/1520-0493(1999)127%3C2204:FASOCR%3E2.0.CO;2).

COAWST's SCRIP-file generator never reads the staggered arrays. Lib/SCRIP_COAWST/read_wrf.f reads only XLONG/XLAT — mass points — and create_psimesh averages four surrounding mass points into psi corners, yielding (nx-1, ny-1), i.e. interior corners only. create_fullgrid.f then reaches the boundary by linear extrapolation:

!extrapolate East boundary
x_full_psi(iend,j) = x_psi(iend-2,j-1) + (x_psi(iend-2,j-1) - x_psi(iend-3,j-1))

one ring out on all four sides, then the four corners from the already-extrapolated edges.

So the comparison is:

	corner source	boundary ring
COAWST	average of 4 mass centres	extrapolated, first order
ERF, with this fix	WRF's actual face coordinates	read from the file

ERF's approach is the better-conditioned one provided the column survives the read — WRF genuinely stores the face positions, so the boundary ring is data rather than an approximation. This change is what stops that advantage being discarded. A COAWST-style extrapolated ring does exist in the driver behind driver.remap_extend_corner_ring, defaulting off, so nothing was masking the defect.

MOM6 offers the closest structural precedent for handling the same problem: extrapolate_metric extends the supergrid latitude (var(i,j) = 2.0*var(i,j+1) - var(i,j+2)) under a sentinel guard, applied before subsampling to geoLatT/Bu/Cu/Cv, with the Coriolis parameter recomputed from the result over the full data domain.

sinPhi_m / cosPhi_m are deliberately left on the old clamp

lat_m is dual-purpose, and its two consumers want different values in that final row. The remap wants the true northern edge latitude; the Coriolis factors want the latitude at the point they are differencing to. sinPhi_m and cosPhi_m therefore keep the cell-domain clamp, and this change is bit-identical for every existing ERF consumer.

That is not merely conservative — it is the position the surrounding ecosystem takes:

• ERF differences f into the outermost face. ERF_MakeMomSources.cpp:377-391 reads 0.5*(sphi_arr(i,j,0) + sphi_arr(i-1,j,0)) for x-momentum, 0.5*(sphi_arr(i,j,0) + sphi_arr(i,j-1,0)) for y, and cphi_arr(i,j,0) directly for z — over mfi.nodaltilebox(0/1/2), which is not restricted at domain boundaries.
• WRF, by contrast, never computes the boundary Coriolis at all. SUBROUTINE coriolis sets i_start = MAX(ids+1,its) and i_end = MIN(ide-1,ite) under specified/nested, with the j mirror, excluding both outermost faces. The same guard appears in perturbation_coriolis and curvature.
• WRF implemented the degraded one-sided stencil and removed it. Still present as a comment: ! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD), over ru_tend(its,k,j) = … *0.5*(f(its,j)+f(its,j)) — the duplicating form, tried and discarded.
• ROMS keeps the two semantics apart explicitly. Same halo routine, deliberately different: latr is exchanged with periodicity off (mp_exchange2d(..., .FALSE., .FALSE., ...)) while f is exchanged with EWperiodic(ng), NSperiodic(ng) on — because wrapping a longitude jumps it 360° while wrapping a Coriolis parameter is correct.

The invariant across WRF, WPS, ROMS, MPAS, ICON, MOM6 and NEMO is that no model lets an approximated latitude reach a Coriolis stencil unguarded: either it never differences f there, or it guards the one place it does. ERF currently does neither, so feeding it a corrected-but-extrapolated ghost would move it outside that invariant. Decoupling the arrays is the right long-term fix — ROMS and MOM6 both show what that looks like — but it is a larger change than this one and does not belong here.

Consequence to be aware of: sin_arr(i,ny) != sin(lat_m(i,ny)) in that single ghost row, flagged in a comment at the site. Note also that solverChoice.variable_coriolis defaults to false, so the ghost-differencing path requires an explicit erf.variable_coriolis in the deck.

Verification, and its limits

• Premise, against a real wrfinput_d01: dimensions and the two coordinate values quoted above.
• Mechanism and fix, at runtime: a standalone AMReX program reproducing the exact convert → resize → & → copy → clamped growntilebox chain at the real grid dimensions, run before and after. Before: lon_m(117,50) = -85.38144, the duplicate. After: -85.0964. Max |after − file| over the entire valid region is 0, so valid-region values are bit-identical and the change is confined to the ghost row/column.
• Compilation: full -c object compile of the changed translation unit at -O3 with MPI and NetCDF, zero warnings; and again under a coupled driver's full define set.
• Not done: ERF was never run end-to-end on a wrfinput deck. A full erf_srclib build was out of reach on the machine available. The fix is verified by compile plus a faithful runtime reproduction of the code path, not by a model run. Reviewers wanting an end-to-end check should exercise a real WRF-input case and confirm lon_m(nx,j) != lon_m(nx-1,j).

Because the change is provably confined to the boundary ghosts, a bit-for-bit comparison of an existing WRF-input regression case should be unchanged; any interior difference would indicate the edge-locality argument is wrong and is the single most valuable thing to check.

Interaction with open pull requests

Checked against every open PR touching either file this change edits
(ERF_InitFromWRFInput.cpp, ERF.H), for textual conflict and for
semantic interaction:

PR	title	shared files	assessment
#3800	Base State Fix WRFInput	ERF_InitFromWRFInput.cpp	Hunks begin at :1688, entirely below this change. Its subject is base-state reconstruction (p_hse/pi_hse/th_hse/qv_hse), stopping the per-(i,j,k) copy of WRF's PB/ALB in favour of evaluating WRF's analytic function at ERF's z_phys_cc. It never reads the coordinate arrays. No interaction
#3884	WRF hydrometeor lateral-boundary advection and nudging	ERF_InitFromWRFInput.cpp	Nearest hunks :449 and :620 straddle :486 without touching it; net shift +5 lines before :486, −2 before :863. Semantically disjoint by construction: it makes no edit to subdomain_to_fill_typed, var_fab.resize, the Box intersection copy, or ixType handling, and all of its init_from_wrfinput logic is name-gated to QCLOUD/QICE/QRAIN/QSNOW/QGRAUP, while this change is name-gated to XLONG_U/XLAT_V. Neither can fire on the other's variables. It replaces the NC_names moisture block at indices 24+, leaving XLAT_V // 22 and XLONG_U // 23 as unchanged context, so the coordinate variables keep their positional indices
#3860	Extend SurfaceLayer for all boundary sides	ERF_InitFromWRFInput.cpp, ERF.H	One line at :1595, converting m_SurfaceLayer to an orientation-indexed container — below all hunks here, so it does not even shift them. Its ERF.H hunks are at :1364 and :1550, over a thousand lines from the :158 doc comment here. No interaction with the coordinate arrays

No conflicts expected in any merge order, and the name-gating means the
separation survives future line movement rather than depending on current
offsets.

Worth noting for reviewers running a realistic case: #3884 adds
erf.use_wrf_bdy_qc_qi and erf.bdy_moist_nudge_type = 3. A fully configured
real WRF case is also the configuration most likely to exercise the corner mesh
this change repairs, since it is what puts genuine spatial structure into the
coupled surface fields.

Follow-ups, not addressed here

• InitFromMetgrid fills lon_m/lat_m from XLONG_M/XLAT_M — genuine mass points, not staggered. The corner-mesh consumer's "already nodal in i" assumption is therefore simply wrong for metgrid-initialised runs, independently of this bug. Documented at the accessor rather than changed; probably deserves its own issue.
• i_hi/j_hi derive from boxes_at_level[lev][0] while var_fab is sized from boxes_at_level[lev][idx]. Identical for a single-box level 0; potentially wrong when num_boxes_at_level > 1. Pre-existing and untouched.
• Decoupling the coupling-side corner coordinates from the dynamics-side latitude arrays, per the ROMS/MOM6 pattern above.

_{Provenance: the ERF, COAWST and driver citations above were read directly in-tree. The WRF, ROMS, MOM6, MPAS, ICON and NEMO comparisons come from a source survey of those repositories and are cited by file and symbol so a reviewer can check them; they were not independently re-derived here.}