Viscous IB force integration reads fd_coeff_* outside the interior #1856

A viscous immersed-boundary case whose body sits within fd_number cells of a domain boundary reads fd_coeff_x/y/z outside the range those coefficients were computed for:

At line 1143 of file src/simulation/m_viscous.fpp
Fortran runtime error: Index '-2' of dimension 2 of array 'fd_coeff_y' below lower bound of 0

s_compute_ib_forces samples the viscous stress fd_number cells out from each interior cell (src/simulation/m_ibm.fpp, master line 1059):

do l = -fd_number, fd_number
    call s_compute_viscous_stress_tensor(viscous_stress, q_prim_vf, dynamic_viscosity, i, j + l, k)

s_compute_viscous_stress_tensor then uses that sample index to look up the coefficient itself:

velocity_gradient_tensor(l, 2) = velocity_gradient_tensor(l, 2) + fd_coeff_y(r, &
                         & j)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j + r, k)

but s_initialize_derived_variables_module allocates them for interior cells only (src/simulation/m_derived_variables.fpp:43-48):

@:ALLOCATE(fd_coeff_x(-fd_number:fd_number, 0:m))
@:ALLOCATE(fd_coeff_y(-fd_number:fd_number, 0:n))
@:ALLOCATE(fd_coeff_z(-fd_number:fd_number, 0:p))

So j + l reaches -fd_number, the lookup is out of bounds, and only the reldebug lanes say so. Every other configuration reads whatever is in front of the array and produces a force that looks plausible.

This is pre-existing on master, not introduced by the PR that surfaced it: #1847 adds a 3D_ibm_flapping_plate example and touches neither m_viscous.fpp nor the force path. It is the first viscous IB case in the suite with a body near a domain boundary, which is why nothing caught this before.

Same class as the finite-difference bug in #1792 (coefficients requested outside the range they were computed for), and worth checking whether any other fd_coeff_* consumer indexes at a sample location rather than a cell location.

Fix

Read the nearest interior cell's coefficients: identical on a uniform grid, and on a stretched grid it is the stencil the boundary cell itself uses.

i_fd = min(max(i, 0), m)
j_fd = min(max(j, 0), n)
k_fd = min(max(k, 0), p)

then index fd_coeff_x(r, i_fd), fd_coeff_y(r, j_fd), fd_coeff_z(r, k_fd).

Checked with gfortran, ./mfc.sh test --no-gpu --debug: without the change 86A2C59B (3D -> Example -> ibm_flapping_plate) reproduces the error above; with it that case and both Prescribed Kinematics cases pass. Going up on #1847 since that is where it blocks CI.