Fix an OpenMP data race in flow_dep_bdy called from moist_tile_loop_2 in solve_em

[pp1457]

Fix an OpenMP data race in flow_dep_bdy called from moist_tile_loop_2 in solve_em

TYPE: bug fix

KEYWORDS: OpenMP, data race, reproducibility, lateral boundary conditions, flow_dep_bdy, tiles, solve_em, thread safety

SOURCE: Liao Yun Yang (National Taiwan University)

DESCRIPTION OF CHANGES:

Problem:
When WRF is built with OpenMP (smpar or dm+sm), flow_dep_bdy is called from
inside the moist_tile_loop_2 !$OMP PARALLEL DO in dyn_em/solve_em.F. Its
writes are correctly bounded by the tile (its:ite, jts:jte), but it reads
from a fixed global row/column that a different tile may be writing
concurrently in the same parallel loop.

In flow_dep_bdy (share/module_bc.F, north-edge branch; the other three edges
are analogous):

      IF (jbe - jtf .lt. spec_zone) THEN
        DO j = max(jts,jbe-spec_zone+1), jtf          ! write index: tile-bounded, correct
          ...
            DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
              i_inner = max(i,ibs+spec_zone)
              i_inner = min(i_inner,ibe-spec_zone)
              IF(v(i,k,j+1) .gt. 0.)THEN
                field(i,k,j) = field(i_inner,k,jbe-spec_zone)   ! <-- READ of a FIXED GLOBAL ROW

jbe = jde-1 is the global domain edge, so jbe-spec_zone is a fixed global row
index unrelated to this tile's jts:jte; i_inner is likewise clamped to global
bounds. A tile in the boundary zone therefore reads a row that another tile owns
and is concurrently updating in rk_update_scalar. Both happen inside the same
PARALLEL DO, with no synchronisation between them.

The consequence is that WRF produces different results on every run from
identical inputs. With the default contiguous static schedule each thread
receives one contiguous block of tiles, so the boundary tile and the tile owning
jbe-spec_zone are usually the same thread and the accesses serialise — which
is why this is normally invisible. Any schedule that interleaves tiles across
threads puts them on different threads and the race fires. Because it is a timing
window it is intermittent, so a clean run does not demonstrate absence.

Solution:
Hoist the call out of moist_tile_loop_2 into its own tile loop,
moist_tile_loop_3. The implicit barrier at the end of moist_tile_loop_2 then
guarantees every tile's write has landed before any tile reads the shared row. One
barrier suffices: within the new loop the rows written are
jbe-spec_zone+1 .. jbe while the row read is jbe-spec_zone, so reads and
writes inside the new loop do not overlap.

The boundary block was already the last statement in the tile-loop body, so the
hoist moves it verbatim — no call arguments, conditions or statements are altered,
only the placement of the barrier. The loop-invariant guards
(config_flags%specified, have_bcs_moist, im .ne. P_QV) stay outside the new
loop so no parallel region is entered when the boundary update is not needed.

Scope of this PR. This is the one call site for which we have end-to-end
before/after evidence, submitted separately per the discussion in #2377. Seven
further call sites (tke_2, chem, tracer x2, scalar x3) have the identical defect
and are addressed in a companion PR. A ninth, zero_grad_bdy, is discussed in
#2377 and deliberately held back — its enclosing loop is the acoustic
small-timestep loop, so hoisting there would add a parallel region to the
innermost time loop and warrants its own performance measurement.

ISSUE:
Part of #2377

LIST OF MODIFIED FILES:
M dyn_em/solve_em.F

TESTS CONDUCTED:

1. Do the mods fix the problem?

Yes. Characterised on v4.4, where flow_dep_bdy is byte-identical to the current
code and is called from moist_tile_loop_2 with the same structure, so the same
reasoning and the same fix apply.

Configuration: CONUS 2.5 km real-data case, specified = .true.,
have_bcs_moist = .false., built dm+sm with GCC 13.3.0 / OpenMPI 5.0.9 and
parallel HDF5 / netCDF-4, run with several OpenMP threads per MPI rank and with
enough tiles per patch that tiles interleave across threads rather than each
thread receiving a single contiguous block. 8 time steps, restart files compared
with diffwrf. To make the schedule selectable at runtime, SCHEDULE(RUNTIME) was
added to the solve_em tile loops as instrumentation only; with
OMP_SCHEDULE=static that instrumented build is bit-identical to the
uninstrumented one.

comparison	before fix	after fix
OMP_SCHEDULE=static,1 vs static	113 differing fields	0
OMP_SCHEDULE=static,1 vs itself	113 differing fields	0
OMP_SCHEDULE=dynamic vs itself	113 differing fields	0
OMP_SCHEDULE=dynamic vs static	113 differing fields	0
fixed build, default static, vs unmodified build, default static	—	0

The last row is the important one for existing users: the change is inert
under the default contiguous-static schedule, so it does not alter results.

The corruption is spatially localised exactly as the code predicts. After only 4
time steps, 10 points of MU_1 differ between two identical runs, and every one
lies in the four rows against the north domain edge (the domain has 1200 rows) —
precisely the jbe branch quoted above:

MU_1 differing points: 10
  j (south_north) = 1195, 1196, 1197, 1198
  i (west_east)   = 1328 ... 1335

2. Verification of this patch against the current tree

dyn_em/solve_em.F was preprocessed exactly as the build does it
(sed -> cpp -> standard.exe -> cpp) and compiled with
gfortran -fsyntax-only -fopenmp against a built module set: 0 errors.
Diffing the preprocessed output against the unpatched file shows the only
changes are the relocated ENDDO / !$OMP END PARALLEL DO and the new
moist_tile_loop_3; no call argument, condition or statement differs.
!$OMP PARALLEL DO / END PARALLEL DO counts balance (52/52).

3. Limits of what was tested

The race was characterised at a single configuration. How the failure rate varies
with thread count or tile size was not measured.

4. Are the Jenkins tests all passing?

Not run on our side — we do not have access. This change is confined to the
moisture path and is bit-identical under the default schedule, so we would expect
the regression suite to be unchanged, but it has not been verified.

RELEASE NOTE:
Fixed an OpenMP data race in the lateral boundary update for moisture variables in
solve_em. flow_dep_bdy writes within its assigned tile but reads a fixed global
row or column that another tile may be updating concurrently, so OpenMP builds
could produce different results from identical inputs. The boundary call is now
made from its own tile loop, so the implicit barrier separates the tile updates
from the shared read. Results are unchanged under the default OpenMP schedule.

[weiwangncar]

The regression test results:

Test Type              | Expected  | Received |  Failed
= = = = = = = = = = = = = = = = = = = = = = = =  = = = =
Number of Tests        : 23           24
Number of Builds       : 60           57
Number of Simulations  : 158           150        0
Number of Comparisons  : 95           86        0

Failed Simulations are: 
None
Which comparisons are not bit-for-bit: 
None

[pp1457]

The regression suite has come back clean on this, so it's ready for review whenever convenient.

Two things are ready on our side whenever you want them:

1. The companion PR covering the remaining seven call sites is prepared against current master and can be opened on request.
2. On zero_grad_bdy: our inclination is a separate PR with its own timing numbers, since its tile loop is the acoustic small-timestep loop and hoisting there adds a parallel region to the innermost time loop. Happy to fold it into the second PR instead if you'd rather have fewer.

No urgency from us — just flagging that neither is waiting on us.

[weiwangncar]

@pp1457 I'd suggest that you started the other PRs. Once the PRs are created, we will review them. Thanks!

[pp1457]

Force-pushed a reworded commit message. The tree is byte-identical to the previously tested commit — same diff, no code change — so nothing should need re-verifying.