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.

[asalmgren]

Code review: PR #3860 — "Extend SurfaceLayer for all boundary sides"

• URL: Extend SurfaceLayer for all boundary sides #3860
• Author: Coleman Kendrick (ckendrick), branch most_walls_cloudchamber
• Size: 59 files changed, +195869 / -830 (mostly golden files; ~2.2k lines of code diff)
• Reviewed at head 96560ae, across MOSTAverage / MOSTStress / SurfaceLayer /
  Diffusion / ERF driver / IO

High

1. Inverted assert kills the stress-rotation path

Source/BoundaryConditions/ERF_SurfaceLayer.cpp:947

AMREX_ALWAYS_ASSERT_WITH_MESSAGE(dir == 0 && m_face.isLow(),
                                 "Stress rotation only supported for zlo face")

zlo has coordDir() == 2, not 0. Since ERF.cpp only permits rotation on zlo,
any run with erf.use_rotate_surface_flux = true now aborts — in release builds
too, since this is ALWAYS_ASSERT. Should be dir == 2 && m_face.isLow().

2. set_rotated_fields writes into the wrong (and wrongly shaped) slots

Source/BoundaryConditions/ERF_MOSTAverage.cpp:422

After the U/V/W/T/Qv/Qr reordering, m_rot_fields[2] is allocated on the
z-nodal zvel BoxArray (line 209) while [3]/[4] are on the CC BoxArray. The
copy does MultiFab::Copy(*rot_fields[2], *fields[3], …), i.e. cell-centered
Theta into the z-nodal fab — a BoxArray mismatch (assert in debug, UB in
release) — and rot_fields[3] receives Qv, so the θ average is computed from qv
and rot_fields[4] is never refreshed at all. Indices should be
rot_fields[3] ← fields[3] and rot_fields[4] ← fields[4].

3. sm_index mutated inside the MFIter loop

Source/BoundaryConditions/ERF_MOSTAverage.cpp:1268 (also :1358, :1471)

On a high face, if (dir != imf) sm_index += 1; — and the unconditional
sm_index += 1 in the Tv and Umag blocks — executes once per box. The first
matching box moves the sampling plane to bigEnd(dir)+1 (a ghost layer for
cell-centred fields), after which every later box fails
vbx.bigEnd(dir) != sm_index and is skipped. The plane average therefore sums a
single box while ncell_plane still assumes the whole plane.

It is also a write to an enclosing-scope variable inside #pragma omp parallel,
i.e. a data race. sm_index must stay loop-invariant.

4. New zlo clauses break (and can crash) the existing stretched-grid MOST path

Source/Diffusion/ERF_DiffusionSrcForState_S.cpp:122 — repeated in all four
branches and for ybx

else if (SurfLayer_on_zlo && qty_index == RhoTheta_comp) xflux = hfx_x(i,j,0);

Two problems, both hitting configurations that work today:

• No rotate guard (unlike the _T variant). For a plain
  zlo.type = surface_layer run on StretchedDz, the horizontal θ flux in the
  lowest cell layer is replaced by SFS_hfx1 — zero or stale without rotation —
  instead of -rhoAlpha*dθ/dx.
• The RhoQ1_comp sibling reads qfx1_x(i,j,0), but SFS_q1fx1_lev is
  nullptr whenever moisture is on and rotation and the X/Y surface layers are
  off (Source/ERF_MakeNewArrays.cpp:677) → null Array4 dereference.

5. m_SurfaceLayer indexed before it is sized

Source/ERF.cpp:756 and Source/Initialization/ERF_InitFromWRFInput.cpp:1598

m_SurfaceLayer.resize(2*AMREX_SPACEDIM) happens at ERF.cpp:1134, inside
InitData_post. But m_SurfaceLayer[Orientation(z,low)] is passed to
update_sst_tsk from init_from_wrfinput — reached via
InitData_pre → InitFromScratch — and from the restart / wrflowinp block at
line 756. Indexing an empty amrex::Vector is out of bounds, and
update_sst_tsk then evaluates if (SurfLayer) on that garbage reference
(Source/IO/ERF_ReadFromWRFLow.cpp:209) → segfault/UB for WRF-input runs with a
low file. Source/ERF_MakeNewLevel.cpp:545 (MakeNewLevelFromCoarse) has the
same exposure. Resize the vector in the ERF constructor.

Medium

6. qr aliased to qv in compute_eb_averages

Source/BoundaryConditions/ERF_MOSTAverage.cpp:2293

(fields[5]) ? fields[4]->const_array(mfi) — the null test was updated to the
new Qr slot but the array still comes from Qv, so
vfac = 1 + epsv*qv - qr collapses to 1 + (epsv-1)*qv. The
virtual-potential-temperature average is wrong for EB + moisture.

7. New Tau* plot variables dereference null Array4s

Source/IO/ERF_Plotfile.cpp:390

tau12 / tau21 / tau31 / tau32 are captured as possibly-empty Array4s and
then used unconditionally in the three ParallelFors.
Tau[lev][tau21/tau31/tau32] is nullptr in the common no-terrain /
no-X-Y-surface-layer case (Source/ERF_MakeNewArrays.cpp:610), and all nine are
null when l_use_diff is false. So erf.plot_vars_1 = … Tau13 crashes instead
of plotting.

8. Same for hfx3, q1fx3, q2fx3

Source/IO/ERF_Plotfile.cpp:437

These three are dereferenced without the (ptr) ? … : 0.0 guard used for
hfx1 / hfx2 / q1fx1 / q1fx2. SFS_q1fx3_lev / SFS_q2fx3_lev are null
for any dry run, so requesting any variable in this group in a dry simulation
dereferences null.

9. Off-by-one on the xhi/yhi tangential stresses

Source/BoundaryConditions/ERF_SurfaceLayer.cpp:873 and :937

For dir == 0, bxx has already been moved to the high x-node
(bxx.setSmall(dir, bxx.bigEnd(dir)), line 826), yet out_ind = k+1 shifts the
z index instead. t31 on the xhi face is therefore written one z-node above
the location used on xlo (out_ind = k). Same pattern at line 937 for t32 on
yhi. The sibling branches get this right by applying the +1 to the face-normal
index (j+1 at line 877, i+1 at line 933).

10. X-face tangential stresses appear swapped

Source/BoundaryConditions/ERF_MOSTStress.H:2255 and :2320
(moeng_flux::compute_u_flux / compute_v_flux, dir == 0)

For dir == 0:

• compute_u_flux builds the stress from velx = v with umean = vm_arr — the
  v-momentum stress — but the caller stores it into t31/t13, the
  w-momentum flux.
• compute_v_flux builds it from vely = w with vmean = wm_arr — w-momentum —
  and the caller stores it into t21/t12, the v-momentum flux.

The Y-face case is self-consistent (u→t12, w→t32), so the X face is the odd
one out: the v and w equations receive each other's wall stress.

11. High-face zref / k-index in set_k_indices_N

Source/BoundaryConditions/ERF_MOSTAverage.cpp:625

lk = bigEnd(dir) - floor((zhi - zref)/dz - 0.5) mirrors an index that is
already correct for a wall-distance zref. With the documented
erf.zhi.most.zref = 0.03125 on the 32-cell/unit-height setup this yields
lk = 1 — the zhi surface layer samples the mean state one cell above the
bottom wall.

Line 627 also stores m_zref = (lk+0.5)*dz, an absolute height, while MOST uses
zref as the distance to the wall in log(zref/z0) and in Rib. Even the
default (zref = ProbHi - dz/2, line 593) then gives zref ≈ 0.98 instead of
dz/2.

12. set_k_indices_T zhi branch pins the k index to a ghost plane

Source/BoundaryConditions/ERF_MOSTAverage.cpp:677 region (~line 852)

In the high-z terrain branch the search loop sets k_arr(i,j,k) = kmax
(= Domain().bigEnd(2)+1) rather than the found lk — the intended line is
present but commented out (//k_arr(i,j,k) = lk;). So compute_*_averages
samples the field in the first ghost cell above the domain regardless of the
requested zref, and zref_arr comes out as dz/2. The low-face branch
correctly uses lk.

Low

13. Sixteen new user-facing plotfile variable names are undocumented

Source/ERF.H (derived_names), Source/IO/ERF_PlotfileSelection.H

Tau11…Tau33, hfx1/2/3, q1fx1/2/3, and q2fx3 are added to the selectable
plot names, but only Docs/sphinx_doc/SurfaceLayer.rst was touched. Per the
project rule these belong in
Docs/sphinx_doc/plotfiles/Plotfile3DReference.rst, including the restriction
that they are only meaningful when the corresponding Tau/SFS MultiFabs are
allocated — which is exactly the condition that findings 7 and 8 turn into
crashes.

The new erf.<face>.most.* / erf.<face>.surface_layer.* input prefixes are
documented.

Lower-confidence observations (not counted as findings)

• compute_fluxes now iterates gtbx over ~4 ghost k-planes for a zlo face:
  redundant MOST iterations, and iterate_flux reads zref/z0 ghost values —
  a NaN there would trip amrex.fpe_trap_invalid.
• The switch from geom.periodicity() to Periodicity(TheDimensionVector(dir))
  in the average FillBoundary calls drops periodic images that the subsequent
  clamp block appears to overwrite anyway.
• const int dir at Source/IO/ERF_Checkpoint.cpp:354 is unused (harmless;
  -Werror is off in CI).

[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.

[asalmgren]

Plane averages are silently zero on every high face.

For any field that is cell-centered in dir — the two tangential velocities (imf < 3 && imf != dir), Theta (imf == 3) and Qv (imf == 4) — this makes sm_index = Domain().bigEnd(dir) + 1. No box's vbx.bigEnd(dir) can ever equal that, so every box hits the continue at line 1311 and the average is left at 0.

The +1 is only correct for the face-normal nodal velocity, which the preceding if (imf < 3 && imf == dir) already handles. Compare compute_region_averages (line ~1749), which computes the same quantity with a plain bigEnd(dir) and no extra increment.

Effect: run the new ABL_MOST_Cloudchamber input with the default most.average_policy = 0 and the xhi/yhi/zhi faces get tm = 0 and umm = WSMIN in surface_temp::iterate_flux, producing garbage u*/t* and wall fluxes. The new test uses average_policy = 1 everywhere, so CI does not cover this path at all.

[asalmgren]

Same class of bug as line 1290. fields[4] (Qv) is cell-centered — the comment three lines above literally says "This is cell-centered so we don't need to worry about double-counting" — so pbx.bigEnd(dir) maxes out at Domain().bigEnd(dir) and can never equal bigEnd(dir) + 1. The virtual-potential-temperature average (iavg = 5) is therefore never accumulated on any high face.

The corresponding region-average block at line ~1749 uses bigEnd(dir) with no +1, which confirms the intended value.

[asalmgren]

Same bug again, for the tangential velocity magnitude. The MFIter is over fields[imf_cc] with imf_cc = 3 (cell-centered Theta), so pbx.bigEnd(dir) != bigEnd(dir) + 1 always and iavg = m_navg-3 .. m_navg-1 (Umag, Umag_xz, Umag_yz) stay zero on high faces.

Region-average counterpart at line ~1870 uses bigEnd(dir).

[asalmgren]

Three problems in the high-face branch.

(a) m_zref is not snapped to the selected cell centre. The low-face branch above stores (lk + 0.5)*m_dz + m_zlo; here the raw user input is stored. With the cloudchamber setup (zref = 0.03125, dz = 0.03125) both faces select the wall-adjacent cell, whose centre is dz/2 = 0.015625 from the wall — but xlo stores zref = 0.015625 and xhi stores 0.03125. log(zref/z0) and Rib then differ, so a geometrically symmetric chamber produces asymmetric fluxes on opposite walls.

(b) The formula disagrees with the asserts and with the default. lk = bigEnd(dir) - floor(zref_tmp/m_dz - 0.5) treats zref_tmp as a distance from the wall, but the two asserts just above (and the default computed at line ~613, ProbHi(dir) - 0.5*m_dz) treat it as an absolute coordinate. So when the user does not set most.zref, the high-face default gives

zref_tmp = ProbHi - dz/2
lk = bigEnd - floor(N - 0.5 - 0.5) = bigEnd - (N-1) = 0

i.e. the xhi/yhi/zhi face samples the cell adjacent to the low wall. The new test only avoids this because it sets most.zref explicitly on all six faces.

(c) Mixed level indexing. m_geom[0].Domain().bigEnd(dir) is combined with m_dz from m_geom[lev], so on a refined level the returned index refers to the level-0 domain.

[asalmgren]

z_lo = face + lk*dx and z_hi = face + (lk+1)*dx, so myhalf*(z_hi + z_lo) = face + (lk+0.5)*dx and this line evaluates to 2*prob_lo + (lk+0.5)*dx.

The high-face branch below (line ~815) does face - myhalf*(z_hi+z_lo), which correctly reduces to the wall distance (lk+0.5)*dx. This line should subtract, not add — it is only correct today because geometry.prob_lo == 0 in the new test.

Separately, k_arr(i,j,k) = k with //k_arr(i,j,k) = lk; commented out just above looks unintended — the whole point of the search loop is to find lk.

[asalmgren]

Disabling this parameter re-introduces double-counting of the surface scalar flux under native SHOC.

The flux override is now driven by SurfLayer_zlo = (SurfLayer[2] != nullptr) (ERF_SetupDiff.H:10) rather than by this flag. But the callers still compute and pass a meaningful value:

• ERF_SlowRhsPre.cpp:267 sets l_apply_surface_layer_fluxes_in_diffusion = false when native_shoc_lev->owns_scalar_surface_fluxes(), precisely so diffusion does not re-apply the MOST heat/moisture flux at k == dom_lo.z;
• ERF_SlowRhsPost.cpp:155-165 does the same for EAMxx SHOC and for native SHOC without uses_host_diffusion().

With the parameter commented out, that decision is silently discarded and zflux(i,j,dom_lo.z) = hfx_z(...) is applied anyway, double-counting the surface scalar flux for surface_layer + native-SHOC runs.

_S.cpp:78 has the identical change. _T.cpp is inconsistent rather than uniform: it still uses use_SurfLayer at line 873, but its SurfLayer_on_zlo checks no longer consult it — and it now declares an unused parameter, which will warn under -Wunused-parameter.

Suggested fix: keep the flag and gate on use_SurfLayer && SurfLayer_zlo.

[asalmgren]

Behavior change on the existing zlo path — please confirm it's intended.

This was mdata.kappa * (qvm_arr(i,j,k) - q_surf_arr(i,j,k)) / (C - psi_h). The new qv_s falls back to qvm_arr whenever q_surf_arr <= 0 (lines 1341-1349), and default_land_surf_moist is 0 (ERF_SurfaceLayer.H:1278, applied at :600).

So for any existing zlo moist run that does not set most.surf_moist, q_star changes from kappa*qvm/(C - psi_h) to exactly 0 — the surface moisture flux is silently removed. The fallback itself looks deliberate (the comment says "there will be no flux"), and it's arguably more physical than the old behavior, but it is a change to previously-working zlo results rather than new-face-only behavior. Worth calling out in the PR description, and worth checking whether any existing moist MOST regression baseline should have moved.

[asalmgren]

ReadCheckpointFileSurfaceLayer() itself loops over all orientations (ERF_Checkpoint.cpp:1631,1638), but this call sits inside the per-orientation construction loop. On restart with N active faces it runs N times and re-reads every already-constructed face's Ustar/Tstar/MOSTAvg* MultiFabs — O(N²) VisMF reads, and the first invocation happens before the remaining faces exist.

Functionally idempotent, but it should move to after the orientation loop closes.

[asalmgren]

Unguarded dereference. The Tau12/Tau13/... paths below carefully fall back to Array4<Real>{} when the pointer is null, but Tau11 (and the Tau22/Tau33 copies) dereference Tau[lev][...] directly.

When l_use_diff == false, all nine entries are set to nullptr (ERF_MakeNewArrays.cpp:714), so erf.plot_vars_1 = ... Tau11 in a run without diffusion is a null dereference. The new doc rows in Plotfile3DReference.rst say "only available when diffusion is enabled", but nothing enforces or diagnoses it — either guard the copy or Abort with a clear message during plot-variable validation.

[asalmgren]

This else (no-diffusion) branch nulls SFS_hfx1/2/3, SFS_q1fx1/2/3 and SFS_diss, but not SFS_q2fx3_lev[lev] — even though the moisture-off branch just above (line 710) was extended to do exactly that.

On a RemakeLevel that turns diffusion off, SFS_q2fx3_lev[lev] can therefore retain a MultiFab built on the previous BoxArray. The new plotfile block at ERF_Plotfile.cpp:1878 reads SFS_q2fx3_lev[lev]->array(mfi) guarded only by a null check, so a stale pointer would be indexed with an mfi from a different BoxArray.