Results depend on the order HRUs and channels are simulated (three conditionally-written SAVEd locals)

[rafiei-vahid]

Three variables in the engine are read on execution paths where they were never written during that
call. Because each is declared with an initializer — and so carries the SAVE attribute under
F2018 8.5.16 — it retains the value from the previous call, and these routines are called once per
HRU or per channel. The result is that model output depends on the order in which spatial objects
are simulated.

This is not a threading issue. Everything below is single-threaded, and the effect is present in the
released code as distributed.

Demonstration

We built main at cb442f7 unmodified, then built a second binary differing only by the two
one-line fixes proposed below, and ran both with OMP_NUM_THREADS=1 on an 11,281-HRU model for two
simulated years with in-stream water quality active:

values compared:  1,245,473,796
values differing:        11,695

  cbod_out / cbod_in / cbod_stor    ~3,500 channel-days each   <- enratio
  chla_out / chla_in / chla_stor      ~600 channel-days each   <- gra
  dox_out  / dox_in  / dox_stor          58 channel-days       <- via cbodu

Adding two lines that can do nothing except remove a carry-over between objects changes the answer.

We should be straightforward about magnitude: chlorophyll-a moves by about 0.2 %, and although CBOD
moves by more than 100 % on those channel-days, CBOD is near zero in this model so the relative
figure overstates its physical importance. The concern is not the size of the change. It is that the
value depends on which object happened to be simulated immediately before, which is not a modelling
choice anyone made and which a user has no way to detect.

1. ch_watqual4.f90 — algal growth rate gra

 22:  real :: gra = 0.
...
186:  if (algcon < 5000.) then
187:    select case (ch_nut(jnut)%igropt)
188:    case (1) ... 191: case (2) ... 194: case (3) ...
201:    end select
202:  end if
205:  factk = Theta(gra,thgra,wtmp) - Theta(ch_nut(jnut)%rhoq, thrho, wtmp)

gra is read unconditionally at 205, but written only inside the algcon < 5000. guard and inside a
select case with no case default. Two reachable paths leave it unwritten:

• algcon >= 5000., which the code reaches by its own construction: line 235 clamps
  algcon_out to 5000 and line 238 exports the corresponding ht3%chla, so the next channel
  downstream recomputes algcon = 5000. exactly at line 165 and the strict < test fails. With
  ai0 = 50, algcon = 20 · chla, so low-flow channels reach the cap routinely.
• igropt outside {1,2,3}. ch_read_nut.f90 clamps only igropt <= 0; values ≥ 4 pass through
  and never write gra on any call.

Affects ht2%chla (via factk → alg_m → algcon_out) and ht2%solp (via zz, line 333).

Suggested fix — one line, consistent with what the cap means, since at saturation there is no
further algal growth:

gra = 0.
if (algcon < 5000.) then

2. varinit.f90 — enrichment ratio enratio is shadowed, so its reset never takes effect

 58:  use hru_module, only : hhqday, ihru, albday, ... , vpd, fixn    ! enratio NOT imported
 71:  real :: enratio = 0.      !none  |enrichment ratio calculated for day in HRU
 92:    enratio = 0.            ! zeroes the LOCAL; the module variable is untouched

varinit resets the HRU's daily state, and its own header comment lists enratio among the
variables it initializes — but enratio is declared as a local, shadowing hru_module's
variable of the same name. The reset has therefore never reached the variable it was written for.

The module variable is written only by pest_enrsb, which hru_control calls only when
surfq(j) > 0. .and. qp_cms > 1.e-6 .and. precip_eff > 0.. On every other HRU-day it retains the
previous HRU's enrichment ratio. Its consumer is swr_subwq.f90:78,
org_c = (soil1(j)%cbn(1)/100.) * enratio * sedyld(j) * 1000. → cbodu(j) → doxq(j), which reach
the channel through the hydrograph. Note the recent carbon rework (dfce092) moved the hsc_d path
behind cswat == 2, which makes this the default branch rather than the alternative.

Instrumenting it on the model above, enratio differed on 28,001 HRU-days between two runs that
visited HRUs in different orders, while sedyld, qdr, surfq, qp_cms, precip_eff and soil
carbon were all bit-identical, and the guard deciding whether enratio is written was false in both.

Suggested fix: add enratio to the use hru_module, only : list and delete the local
declaration.

etday, crk, over_flow and sedprev are shadowed in varinit the same way. etday is written
unconditionally in hru_control and the other three are unused there, so they appear harmless — but
they are worth removing for the same reason.

3. sd_channel_control3.f90 — storage coefficient scoef

 41:  real :: scoef = 0.
265:  if (rttime > time%dtm / 60.) then       ! travel time > routing time step
267:    scoef = 24. / (ch_rcurv(jrch)%in2%ttime + ch_rcurv(jrch)%out1%ttime + 24.)
...
283:  else
285:    hdsep2%flo_surq = scoef * hdsep1%flo_surq     ! never assigned on this path
...
348:    hcs2%salt(isalt) = scoef * hcs3%salt(isalt)   ! read outside the if entirely

scoef is assigned only in the then branch. The else branch — taken whenever travel time is
shorter than the routing step, i.e. under 24 h at a daily step, which is the normal case for most
reaches — multiplies by whatever the last long-travel-time channel left behind. The salt and
constituent lines from 348 are read outside the conditional entirely. In a model where no channel
exceeds 24 h travel time, scoef is never assigned at all.

Affects ob(icmd)%hdsep%flo_* (hydrograph separation) and channel salt/constituent outflow.

We are not proposing a fix for this one. The comment on the else branch reads "route all stored
and frac of incoming", which suggests a different coefficient was intended rather than the same one,
and we would rather you decide what it should be than guess. We also have no differential reproducer
for it, since it fails identically on every path we can construct — the report rests on the control
flow above.

The common mechanism, and auditing the rest

All three share a cause that is a property of Fortran rather than of the hydrology: a local declared
with an initializer acquires the SAVE attribute and persists across calls. No compiler option
changes this — -auto, -recursive and -qopenmp all exclude SAVEd variables by definition.

Commit 39fabde (2024-08-08, "Initialized varables with python script...") added = 0. to
declarations across the code base. That is ordinarily good hygiene, and in most languages it would
be; in Fortran it made those locals persistent. Initialized declarations in src/*.f90 went from
2,353 to 8,024 in that commit, and every conditionally-written one among them became a potential
carry-over between objects. Both defects we can demonstrate belong to that population.

For scale, we scanned other community models for the same pattern: SWAT 2012 has 66 initialized
locals inside procedures, SUMMA 14, mizuRoute 8, against roughly 3,400 in current SWAT+.

We are happy to share the scanner we used, and to open PRs for items 1 and 2 whenever you would like
them — they are one line each, and we have kept them separate from any other work.

---

For completeness on how these surfaced: we were building a multi-threaded version of the engine and
required it to reproduce single-threaded output bit for bit. That requirement is what made the
carry-overs visible. The defects themselves are independent of it, as the single-threaded
demonstration above shows.