Enable refinement with STF/Sullivan terrain smoothing

[adamwise95]

Summary

This PR enables adaptive mesh refinement (AMR) with terrain smoothing methods STF (terrain_smoothing=1)
and Sullivan (terrain_smoothing=2). Previously, ERF would abort when attempting to use multiple refinement
levels with non-BTF terrain smoothing.
Key changes:

• Interpolate terrain coordinates from coarse to fine levels for STF/Sullivan methods
• Add two refinement strategies: "interpolate" (default) and "transform" modes

Implementation Details

Interpolate Mode (default: erf.amr_terrain_refinement = "interpolate")

• Fine levels (lev > 0) use terrain coordinates interpolated from the coarse level
• Maintains consistency between levels
• Simpler and more conservative approach

Transform Mode (erf.amr_terrain_refinement = "transform")

• Fine levels read high-resolution terrain data from file
• Blends with interpolated coordinates using height-dependent decay
• Captures finer terrain features

Modified Files

• Source/ERF_MakeNewArrays.cpp - Add interpolation vs transform logic for fine levels
• Source/Utils/ERF_TerrainMetrics.cpp - Early return for fine levels with STF/Sullivan
• Source/Utils/ERF_PlaneAverage.* - Fix multilevel planar averaging
• Source/BoundaryConditions/ERF_MOSTAverage.cpp - Fix surface layer for late-start levels

Testing

Tested with:

• 2-level AMR with terrain_smoothing=1 (STF)
• Both "interpolate" and "transform" refinement modes
• Various terrain configurations

Related Issues

Resolves the limitation that prevented AMR with STF/Sullivan terrain smoothing.

Additional Bug Fixes

• Fix multilevel planar averaging compatibility
• Fix surface layer initialization for levels starting at t > 0 (for surface cooling)
• Changed abort when terrain heights are taller than the specified height to be a warning.

[asalmgren]

Review: PR 3935 — Enable refinement with STF/Sullivan terrain smoothing

Branch STF_refinement_PR (6 commits). Files changed: Source/ERF_MakeNewArrays.cpp,
Source/Utils/ERF_TerrainMetrics.cpp, Source/Utils/ERF_PlaneAverage.H,
Source/BoundaryConditions/ERF_SurfaceLayer.H, Docs/sphinx_doc/Inputs.rst.

The core issue: what the removed abort was protecting

make_terrain_fitted_coords opens with

z_phys_nd.setDomainBndry(bogus_large_value,0,1,geom);   // ERF_TerrainMetrics.cpp:62

and depends on init_which_terrain_grid to overwrite every domain-external ghost
cell afterwards — the lateral ghosts via the clamped ii/jj loop over growntilebox,
and the bottom layer via z_arr(i,j,-1) = 2*z(i,j,0) - z(i,j,1).

This PR turns that function's Abort("Must use terrain_smoothing = 0 when doing multilevel") into an early return. Consequently, for any lev > 0 with
terrain_smoothing != 0, those ghost cells remain at bogus_large_value (1e18/1e150).
FillBoundary cannot repair them in non-periodic directions, and both make_zcc (reads
z_nd over growntilebox) and get_dzmin_terrain (reads z_nd(i+1,j+1,k) at valid-box
edges) consume them — so the poison propagates into z_phys_cc's ghost layers and into
the dzmin handed to the microphysics model.

The PR compensates in only one place (the k=-1 layer, only in init_zphys), which
leaves two other call paths broken and one of the two compensations dead.

Blocking findings

1. Source/Utils/ERF_TerrainMetrics.cpp:188 — early return skips the ghost repair

As above. The function poisons the ghost cells itself and no longer un-poisons them.

2. Source/ERF_MakeNewArrays.cpp:942 — regrid path has no compensating fill

remake_zphys calls InterpFromCoarseLevel with IntVect(0,0,0) (deliberately not
filling outside-domain ghosts) and relies entirely on the following comment being true:

// This recomputes the fine values using the bottom terrain at the fine resolution,
//    and also fills values of z_phys_nd outside the domain
make_terrain_fitted_coords(lev,geom[lev],*temp_zphys_nd,zlevels_stag[lev],phys_bc_type);

With the early return it is no longer true, and no fix-up was added here. The first
regrid of a 2-level STF run therefore leaves z_phys_nd(i,j,-1) and the lateral
boundary ghosts on the fine level at bogus_large_value — even in the default
"interpolate" mode.

Separately: in "transform" mode the high-resolution surface and blended interior
computed in init_zphys are silently discarded at every regrid (pure coarse
interpolation is restored), so the terrain changes discontinuously in time.

3. Source/TimeIntegration/ERF_TI_utils.H:160 — MovingFittedMesh at lev>0

The abort was the only guard on this path. Each step it copies the new surface into k=0
of z_phys_nd, z_phys_nd_src and z_phys_nd_new, then calls
make_terrain_fitted_coords, which now no-ops. All k>0 nodes keep stale values — for
_src/_new, essentially uninitialized apart from the k=0 slab — and the subsequent
make_J calls turn that into garbage detJ and z_t_rk. The same applies to the
fine-level calls in ERF_InitFromWRFInput.cpp:1277 and ERF_InitFromMetgrid.cpp:338,
which the new docs declare unsupported but which nothing now rejects.

4. Source/ERF_MakeNewArrays.cpp:860 — transform-mode k=-1 fill is dead code

if (k0 == domlo_z) {
    ...
    if (k0 > 0) {              // never true: domlo_z == 0 for ERF domains
        // fill z_arr(i,j,k0-1)
    }
}

The two conditions are mutually exclusive, so z_arr(i,j,-1) is never written. Combined
with the setDomainBndry poisoning, every "transform"-mode run leaves the entire
k=-1 node layer of the fine level at bogus_large_value, which make_zcc copies into
z_phys_cc's lower ghost cells. The interpolate-mode block at line 890 does fill it,
which shows the intent.

5. Source/ERF_MakeNewArrays.cpp:840 — decay profile depends on box decomposition

Real z_patch_top = zlevels_stag[lev][khi];   // khi = bx.bigEnd()[2] of THIS box

Two bottom-touching boxes with different vertical extents — which STF explicitly
tolerates, cf. the all_boxes_touch_bottom re-chop for lev 0 — get different decay
profiles. Nodal boxes share their boundary node plane, so those shared nodes are assigned
two different values in two FABs and FillBoundary picks one arbitrarily: a mismatched
grid at the seam. Even with no seam, changing max_grid_size or the rank count changes
the computed terrain.

6. Source/ERF_MakeNewArrays.cpp:764 — amr_terrain_refinement is unvalidated

Any value other than exactly "interpolate" or "transform" (a typo, "Transform",
"none") leaves both flags false and takes the original path: the fine terrain is read
into k=0, init_which_terrain_grid returns early so no vertical transform is applied,
and neither post-fix block runs. The result is k=0 at fine-resolution terrain with
k=1..N coarse-interpolated — non-monotonic (negative detJ) wherever the fine surface
exceeds the interpolated k=1 height — with no diagnostic at all. Needs an Abort
on unrecognized values, and one for init_type == WRFInput/Metgrid, which the new docs
call unsupported.

Other findings

7. Source/ERF_MakeNewArrays.cpp:890 — interpolate-mode fill reads poisoned ghosts

The extrapolation reads z_arr(i,j,0) and z_arr(i,j,1) over growntilebox, i.e.
including outside-domain cells — exactly the ones setDomainBndry left at
bogus_large_value. So z_arr(i,j,-1) becomes bogus too wherever the fine level touches
a non-periodic lateral boundary. The fix works only for a laterally interior or fully
periodic fine level.

8. Source/ERF_MakeNewArrays.cpp:906 — "Terrain is taller than domain top!" suppressed

The abort is commented out (with the code comment conceding "This needs to be fixed in
the STF algorithm") and replaced by three unconditional amrex::Print lines on every
level-0 init. When zmax > z_top, the BTF/STF map

((z_sfc - z_lev_sfc)*z_top + (z_top - z_sfc)*z) / (z_top - z_lev_sfc)

has negative slope in k, so z_phys_nd decreases with height and detJ/dz go
negative — the run continues and produces NaNs far from the cause. It also interacts with
the new transform blending, whose monotonicity only holds while
delta_terrain < z_patch_top. Prefer a tolerance for the known small STF overshoot over
removing the check; drop the unconditional prints.

9. Source/Utils/ERF_PlaneAverage.H:287 — empty-plane guard defers the failure

Suppressing the division leaves m_line_average at the accumulated sum, i.e. exactly 0,
rather than marking the plane invalid. ERF_MakeSources.cpp:373-374 then computes
dptr_t_plane(k+1) / dptr_r_plane(k+1) for the subsidence gradient; for the topmost real
cell of a fine patch that does not span the domain height, plane k+1 has zero cells, so
this is 0/0 -> NaN (and 0 - 300 -> a huge spurious gradient in mixed cases).
m_ncell_plane = cnt_h[0] likewise becomes 0 for a level whose grids do not reach the
bottom. Empty planes need a nearest-valid fill or an explicit sentinel consumers check.

10. Source/BoundaryConditions/ERF_SurfaceLayer.H:911 — right direction, still unsafe

Switching m_geom.size() to t_surf.size() fixes the out-of-bounds read for the
late-start-level case, but the vector is resized to nlevs (finest_level+1, see
ERF.cpp:1285 / ERF_MakeNewLevel.cpp:851) on the first per-level call, so entries for
levels not yet built are nullptr and t_surf[lev]->setVal dereferences null if
update_surf_temp runs in that window (or after a level is cleared). Add
if (t_surf[lev]); same pattern for the sibling arrays resized alongside it.

11. Source/ERF_MakeNewArrays.cpp:816 — leftovers in the transform block

domhi_z and z_top are computed and never used; nz is used only in the redundant
k < nz test (valid nodal k never exceeds nz-1). Harmless today (CI has -Werror
commented out), but they read as remnants of a different decay formulation and obscure
which "top" the decay is meant to reference — the ambiguity behind finding 5.

12. Source/ERF_MakeNewArrays.cpp:842 — validbox() inside a tiled MFIter

Both new loops use const Box& bx = mfi.validbox() as the ParallelFor range inside a
TilingIfNotGPU() MFIter, so on tiled CPU builds each box's work is repeated once per
tile. Results are unaffected (writes are idempotent, no #pragma omp parallel here, so no
race), but mfi.tilebox() is the surrounding convention.

Non-findings (checked and cleared)

• delete z_phys_interp / z_lev_d going out of scope immediately after the MFIter loop
  is not a GPU async hazard — MFIter's destructor stream-syncs.
• Documentation: erf.amr_terrain_refinement is added to the Inputs.rst parameter table
  with prose for both modes, so the "document new user-facing options" requirement is met.
  The prose claims WRFInput/Metgrid are unsupported, which the code neither enforces nor