Extend SurfaceLayer for all boundary sides

[ckendrick]

This allows the MOST/Surface Layer boundary condition to be specified on all 6 faces instead of only zlo.

Main changes to SurfaceLayer/MOSTAverage classes:

• This adds a orientation parameter to the SurfaceLayer and MOSTAverage classes to determine behavior on different faces.
• MOST configuration options are prefixed with the face name if the SurfaceLayer BC type is specified on more than the zlo face (e.g, erf.most.surf_temp, can now be erf.xlo.most.surf_temp).
• MOST Averages and SurfaceLayer MultiFabs are constructed on planes corresponding to the face rather than at k=0.
• MOSTAverage is extended to compute averages of W and other planar velocities (XZ, YZ)
• MOSTAverage field ptrs are re-ordered from U/V/T/Qv/Qr/W to U/V/W/T/Qv/Qr.
• Local and regional averaging are supported, but not terrain interpolation, time windowing, EB, etc.
• SurfaceLayer writes to Tau and heat/moisture flux arrays for each face. Directional velocity averages and orientation is passed to MOSTStress.
• Only the Surface Temperature and Moeng flux pathway are updated for now to handle different faces.

Main changes to ERF:

• ERF now constructs a SurfaceLayer instance all 6 faces where the BC type is set to surface_layer (otherwise nullptr)
• Existing coupling dependent on reading surface layer values (SHOC, PBL, etc) should default to checking and reading on the zlo face only. Existing behavior on the zlo face should not be modified by these changes.
• X and Y hfx and qfx are read in the diffusion routines (N/S/T versions) from the surface layer on the face (if set).
• Each face writes its variables to checkpoint with a face index prefix ("UStar_0, UStar_1", etc). Existing checkpoints without a face prefix should be supported.

A follow-up PR will be made to adapt this for the new cloudchamber interface and problem configuration.

[ckendrick]

@asalmgren I think the review comments should be addressed now

[asalmgren]

Review: index handling on the new high faces

Thanks for this — the overall structure (per-orientation SurfaceLayer instances, face-prefixed inputs, face-indexed checkpoint entries) reads well, and the zlo path is mostly preserved. I focused on the new non-zlo code paths and found a consistent class of index bug there, plus one regression that reaches the existing zlo path.

The most important structural point: the new ABL_MOST_Cloudchamber test sets most.average_policy = 1 on all six faces and sets most.zref explicitly on all six, so the plane-average path and the default-zref path on high faces get no CI coverage — and both are broken (comments 1-4 below).

High severity

1. ERF_MOSTAverage.cpp:1290, 1381, 1483 — sm_index = bigEnd(dir) + 1 compared against cell-centered boxes means the plane averages for U/V/Theta/Qv, theta_v and Umag are never accumulated on any high face and end up zero.
2. ERF_MOSTAverage.cpp:645 — high-face zref is not snapped to the cell centre (so opposite walls disagree for a symmetric domain); the lk formula treats zref as a wall distance while the asserts and the default treat it as an absolute coordinate, so the no-most.zref default selects lk = 0; and m_geom[0] is mixed with m_geom[lev].
3. ERF_SurfaceLayer.cpp:976 — bxy is not nodal in x, so the xhi wall stress on w-momentum is written to an interior node and never reaches the boundary node.
4. ERF_SurfaceLayer.cpp:956, 1016, 1020 — the is_low_face ? n : n+1 shift is applied to tangential indices for dir == 0/dir == 1, offsetting opposite walls by one cell.
5. ERF_DiffusionSrcForState_N.cpp:76 (and _S.cpp:78, _T.cpp:94) — use_SurfLayer is now ignored, so the native-SHOC "don't re-apply the surface scalar flux" decision made in ERF_SlowRhsPre.cpp:267 is discarded and the flux is double-counted.

Medium
6. ERF_MOSTStress.H:1417 — q_star becomes exactly 0 for existing moist zlo runs that don't set most.surf_moist. Looks intentional, but it changes previously-working zlo behavior.
7. ERF_SurfaceLayer.cpp:395 — compute_fluxes expands the box to the full domain height while the FABs only own their own box's k-range; out-of-bounds if the grid is decomposed in z.

Low: ERF.cpp:1336 (O(N^2) checkpoint reads on restart), ERF_Plotfile.cpp:1820 (unguarded Tau11 deref with diffusion off), ERF_MOSTAverage.cpp:797 (2*prob_lo offset in the low-face terrain zref), ERF_MakeNewArrays.cpp:717 (SFS_q2fx3_lev not nulled).

Suggested gate before merge: switch the test to most.average_policy = 0 (the block already commented out at lines 74-79 of ABL_MOST_Cloudchamber.i) and drop the explicit most.zref lines. That alone exercises items 1 and 2.

Review produced with Claude Code.

[ckendrick]

@asalmgren The new review issues should be addressed now.
The deviations from the existing zlo path are fixed and checks are added to ensure other setups (SHOC/PBL) use MOST only on the zlo surface. The new ABL_MOST_Cloudchamber test setup is intentionally unchanged to mirror a real setup (using regional averaging and explicit zref on each side). However, unit tests were added to check for both planar and regional averaging with default zref so the missing coverage should be picked up in CI now.

[AMLattanzi]

Summary

Things I'd want addressed before merge

1. Plane averaging double-counts on lateral faces under tiling.

compute_plane_averages, first loop (ERF_MOSTAverage.cpp:1213-1228): the face test is on vbx (the grid), then every tile of that grid gets pbx.setSmall(dir,sm_index); pbx.setBig(dir,sm_index). TileNoZ() returns (mfiter_tile_size[0], mfiter_tile_size[1], 1024000), and AMReX's default mfiter_tile_size is (1024000,8,8) — so on CPU there is tiling in y. For a ylo/yhi face, every y-tile of a grid collapses onto the same y-plane and gets summed, while ncell_plane counts the plane once. A grid 32 cells deep in y over-counts by 4×. x-faces are safe only because the default x tile size is huge; set fabarray.mfiter_tile_size and they break too. The Tv block and the Umag block right below both test pbx instead and are correct — that internal inconsistency is the tell. Switching the first loop to the pbx test matches them and fixes it.

The new gtest wouldn't catch this: the domain is (3,5,7)–(5,8,11), one box, 4 cells in y, so no tile split. Bumping that domain past 8 cells in y (or adding a multi-box BoxArray) would cover it, and would also give the first multi-grid coverage of these paths.

2. t31_arr / t32_arr null-deref guard is inverted. ERF_SurfaceLayer.cpp:~690:

cpp
auto t31_arr = (dir == 0 || Tau_lev[TauType::tau31]) ? Tau_lev[TauType::tau31]->array(mfi) : Array4<Real>{};

When dir == 0 the short circuit skips the pointer check entirely, so a null tau31 dereferences. In practice update_diffusive_arrays allocates tau21/31/32 whenever l_Surf_X || l_Surf_Y, so this only fires with l_use_diff == false — but then it's a segfault where an assert would be much kinder. (Tau_lev[...] != nullptr) alone is the right condition; the dir == 0 clause buys nothing.

3. Ghost values on lateral faces are never communicated.

update_fluxes now skips FillBoundary on u_star/t_star/q_star/olen unless coordDir() == 2. But compute_u_flux for dir == 0 reads u_star_arr(ic,j-1,kc) and umm_arr(ic,j-1,kc) at the low y-node of every grid. With more than one grid in y along an x-wall, that's a ghost cell filled only by whatever compute_fluxes happened to compute over its own growntilebox, using average-MF ghosts that (see 4) are themselves not filled tangentially. The comment says the faces are "filled individually on corresponding ranks," which is true of the valid region but not the halo. Worth either restoring FillBoundary for all faces or documenting why the grown-tilebox computation is sufficient.

4. Periodicity(IntVect::TheDimensionVector(dir)) on the average FillBoundary calls (lines 1677, 1800, 1944-1947).

Periodicity takes period lengths, so this asks for period 1 in the collapsed direction and non-periodic in the other two. It does propagate the single valid plane into the z-ghosts, which I assume is the intent — but it also drops the tangential periodic images that geom.periodicity() used to supply. For a zlo run with periodic x/y and average_policy = 1/2, the ghost values feeding um_arr(i-1,...) at the periodic boundary now come from nowhere. No gold file changed, so either the affected configs aren't in CI or the ghosts don't reach the answer — but the reasoning should be in a comment rather than left implicit.

6. Cross-BoxArray MFIter indexing.

compute_fluxes and fill_qsurf_with_qsat now iterate *m_lmask_lev[lev][0] and call ->array(mfi) on MultiFabs built on a different BoxArray (collapsed in dir, not in z). This works only because both derive from the same 3D grids with the same DM and because of the new full-z-column assert — AMReX will not catch a violation, even in debug. The rationale in the comment is good; I'd add an AMREX_ASSERT on boxArray().size() and DM equality so a future decomposition change fails loudly. Relatedly, gtbx.grow(2,3) in fill_qsurf_with_qsat hardcodes a ghost width that must not exceed ng[0] of the 2D MF — worth deriving from nGrowVect().