statsout: SGS statistics defects survive the rewrite (wrong sign, wsgs is a tendency, dead TKE budget, corrupted bottom level)

[mvreeuwijk]

Companion to #306, which fixes these defects on master. Filed separately because the statsout rewrite (PR #275, destined for v3.0) relocated the SGS code verbatim rather than rewriting it, so none of the master patches will apply — or even textually locate their target — on this branch.

Line references are against origin/statsout @ 75a6e09.

Why a separate issue

statsout splits modstatsdump.f90 in two:

• src/out_stats.f90 (new, module stats) — all of tdump/xydump/xytdump/ydump/ytdump, including every SGS flux. Called from program.f90:211.
• src/modstatsdump.f90 (shrunk) — retains only lmintdump and the ltkedump budget.
• src/modstatistics.f90 was deleted.
• modmpi.f90 gained spatial_avg_xy/spatial_avg_y, near-verbatim copies of avexy_ibm/avey_ibm — including the IBM kb fixup, which now exists in two places.

src/out_stats.f90:1375-1394 is character-for-character identical to master's modstatsdump.f90:835-857, down to the trailing ! tg3315 check this. Of the eleven defects confirmed on master, nine survive unchanged, one is partly fixed, and one (metadata) is mostly fixed — and the rewrite introduced two regressions.

1. HIGH — every SGS flux has the opposite sign to the resolved flux in the same file

src/out_stats.f90:1377, :1383, :1464

diffu (src/modsubgrid.f90:707-729) adds +d/dz[emom*(du/dz + dw/dx)] to up. Since du/dt = ... - dtau13/dz, the physical stress is tau13 = -emom*(du/dz + dw/dx). The statistics store the un-negated quantity:

! out_stats.f90:1377
usgs(i,j,k) = emom * ( (um(i,j,k)-um(i,j,k-1))*dzhi(k) &
                     + (wm(i,j,k)-wm(i-1,j,k))*dxhi(i) )

i.e. usgs = +2*nu_t*S13 = -tau13. Same for vsgs (:1383) and for thlsgs/qtsgs/svsgs via stats_compute_sgs (:1464-1465), which reproduces diffc's face flux (modsubgrid.f90:581) without negating it. The emom weights sum to 1 on a stretched grid, so the magnitude is correct — only the sign is wrong.

The resolved fluxes written into the same file with the same units use the physical convention: out_stats.f90:1580 (upwpxytik = <uw> - <u><w>), :1648, :1607.

Impact: upwp + usgs yields <u'w'> - tau13 instead of the total stress. Affects usgs, vsgs, thlsgs, qtsgs, sNsgs across stats_t, stats_xy, stats_xyt, stats_y, stats_yt. Nothing in docs/ states a convention.

Fix: negate at :1377, :1383, :1464; document in the ncinfo long_names and in docs/udales-output-files.md.

2. HIGH (regression) — wsgs is a tendency labelled as a flux, now written to five files instead of one

src/out_stats.f90:1392-1394

wsgs(i,j,k) = ( ekm(i,j,k)  *(wm(i,j,k+1)-wm(i,j,k))  *dzfi(k)   &
               -ekm(i,j,k-1)*(wm(i,j,k)  -wm(i,j,k-1))*dzfi(k-1) ) * 2. &
               * dzhi(k) ! tg3315 check this

Character-for-character the vertical term of diffw. Dimensions (m^2/s)*(m/s)/m/m = m/s^2, not m^2/s^2. The intended tau33 = -2*nu_t*dw/dz is cell-centred; this is d/dz(2*nu_t*dw/dz) at a w-face — wrong quantity and wrong staggering.

This is a regression against master. Master declared wsgs in exactly one ncinfo; statsout declares and writes it in five, including the 3-D per-CPU dump:

File	Location	Grid code
tdump (3-D)	out_stats.f90:723	ttmt
xytdump	:912	mt
xydump	:1004	mt
ytdump	:1078	t0mt
ydump	:1179	t0mt

A quantity carrying its author's own ! tg3315 check this and the wrong physical dimension has been promoted to a first-class output of four more files.

Fix: wsgs(i,j,k) = -2.*ekm(i,j,k)*(wm(i,j,k+1)-wm(i,j,k))*dzfi(k), with cell-centred grid codes (tttt/tt), or drop the variable if tau33 is not wanted.

3. CRITICAL (regression) — the ltkedump TKE budget writes zeros; its compute engine was deleted

src/modstatsdump.f90:188, src/modfields.f90:480-574

modstatsdump.f90:187-188 still calls tkestatsdump, but src/modstatistics.f90 was deleted on this branch (git diff --stat master...HEAD shows -398 lines). That file held genstats and tkestats — the only code that ever assigned uav, tvmx, tsgsmx1/2, strain2av, nusgsav, disssgsav. All inputs are allocated and zeroed at modfields.f90:527-532, :552-553, :566-567, :571-572, under a comment that already reads ! this is currently not usable.

Result: tkedump.<expnr>.nc is written full of exact zeros with no warning.

On master this was recoverable by uncommenting one line (see #352); here the producer must be restored from history or rewritten.

Three defects are latent inside this path and surface the moment it is revived:

• strainav2 is a scalar used across loop nests. Declared real strainav2 (modstatsdump.f90:238), assigned per-gridpoint at :303/:308/:318/:328 inside the nest :243-342, then read at :419 (disssgsfl(i,j,k) = 2.*nusgsav(i,j,k)*strainav2) inside the separate nest :385-548, which never reassigns it. Every point uses the strain of the corner point (ie,je,ke). Fix: real :: strainav2(ib:ie,jb:je,kb:ke).
• Wrong mean field in the z-component of SGS transport. modstatsdump.f90:434-435 pairs tsgsmz1/2 with vav; must be wav(i,j,k)/wav(i,j,kp). The x- and y-blocks (:428-432) correctly use uav and vav, confirming a copy-paste.
• Duplicated and misindexed shear production. :491-492 writes the <v'v'>d<v>/dy term twice. :503-504 uses vpwpav(ip,j,k)/vpwpav(ip,jp,k) where the paired wav differences step in k; should be vpwpav(i,j,kp)/vpwpav(i,jp,kp).

Fix: either restore modstatistics.f90 and fix the three defects, or delete the ltkedump path outright along with the modfields.f90:480-574 allocation block, the ncinfo placeholders at modstatsdump.f90:86-93, and the namelist flag. Given v3.0 is a breaking release, deleting a silently-zero-filled output is probably the better call — worth deciding here rather than on master (#352).

4. MEDIUM — bottom level of every masked profile is corrupted in IBM runs, and the fixup now exists twice

src/modmpi.f90:762-765 (new spatial_avg_xy) and :664-667 (avexy_ibm)

src/modibm.f90:2177 unconditionally zeroes the bottom masks under libm (IIw(:,:,kb) = 0, IIuw(:,:,kb) = 0, IIvw(:,:,kb) = 0), so IIws(kb) = IIuws(kb) = IIvws(kb) = 0. Both averaging routines then hit:

if ((.not. lnan) .and. (IId(kb)==0)) then
  averl(kb) = sum(var(:,:,kb))   ! unmasked sum, includes solid points
  IId(kb)   = IId(ke)            ! divided by the fluid count at the domain TOP
end if

Every SGS call passes lnan=.false. (out_stats.f90:1599-1601, :1610, :1619, :1653-1655, :1664, :1673), so usgs/vsgs/wsgs/thlsgs/qtsgs at k=kb are sum_all / N_fluid(ke) — a meaningless number, not flagged as missing. The same applies to w (:1576, :1629), upwp (:1580), vpwp (:1581), uw (:1585), vw (:1586), ww (:1590), wpthlp (:1607), wthl (:1608), wpqtp (:1616), wqt (:1617).

Note for whoever ports the master fix: patching avexy_ibm alone leaves the entire out_stats.f90 stack broken. Both copies need it.

Fix: pass lnan=.true. for the statistics calls so k=kb reports -999., and confine the fixup to the modthermodynamics caller it was written for (see the tg3315 22.03.19 comment).

5. MEDIUM — SGS statistics ignore the wall model entirely

src/out_stats.f90:1375-1467; solver side src/modwallfunctions.f90 and src/modibm.f90:990/:1033/:1075/:1120

At walls and IBM facets the solver cancels the interior diffusive flux and substitutes the wall-function stress (diffu_corr, diffv_corr, diffw_corr, diffc_corr). The statistics apply the interior formula everywhere.

• Averaged profiles mask those faces out, so the profile has no surface contribution and the stress budget cannot close at the wall.
• 3-D stats_t writes usgst/vsgst/wsgst/thlsgst/qtsgst/svsgst(:,:,kb:ke) with no mask at all, so those points carry interior-formula values unrelated to the stress the solver applied. (On master this is being fixed by masking; see SGS statistics carry no surface/wall-model flux contribution #353 for the wall-flux accumulation itself.)

6. LOW — uninitialised reads at k = ke+kh

src/out_stats.f90:1389, :1463

wsgs (declared kb:ke+kh at :652) and thlsgs/qtsgs/svsgs (:660, :666, :674) are filled only do k=kb,ke but read at ke+kh by stats_compute_tavg (:1499, :1512, :1521, :1536) and by spatial_avg over kb:ke+kh (:1601, :1610, :1619, :1655, :1664, :1673) — which feeds an MPI_ALLREDUCE.

Partial credit: the rewrite did fix master's sv1k oversight (stats_interpolate_k at :1447-1455 correctly runs kb..ke+kh), and all NetCDF writes are truncated to kb:ke, so nothing corrupt reaches a file. But it will trap under -ffpe-trap=invalid in a Debug build and poisons the reduction with NaN.

Fix: extend both loops to ke+kh, or zero the arrays at allocation.

7. LOW — leftover metadata errors

The rewrite's metadata sweep fixed most of the master issues (vsgs to m^2/s^2, sNsgs to g/m^2s, vpwp to m^2/s^2, w*yt to K m/s) — credit where due. Remaining:

• src/modstatsdump.f90:86-93 — placeholder long_names survive verbatim: 'p_bant production or consumption term' (mangled "buoyant"), 'total viscous transport (?)', 'Total turb???'.
• src/out_stats.f90:811 — PSS is a percentage (computed at :1547) but labelled 'gm/s'.
• wsgs units 'm^2/s^2' at :723, :912, :1004, :1078, :1179 are wrong as a consequence of item 2.

8. LOW — O(dt) staleness in the sampled state

program.f90: subgrid (:153) sets ekm/ekh from the substep input state; tstep_integrate (:204) advances to the output state; stats_main (:211) samples um/vm/wm/thlm at the output state with stale ekm/ekh; boundary (:214) runs afterwards. halos (:206) refreshes lateral halos only, so the kb-1/ke+1 ghost planes read at k=kb are one timestep old. Unchanged from master.

Verified correct — please do not re-litigate

I traced these and found no error on statsout:

• Pointwise-then-averaged ordering is correct throughout — the nonlinear closure is evaluated on instantaneous fields (stats_interpolate_and_sgs_*, :1343-1423) before time/space averaging.
• Staggering and stretched-grid metrics are correct for everything except wsgs. usgs sits at the tau13 corner, vsgs at tau23, scalars at the w-face; the emom interpolation weights sum to 1 on a stretched grid; dzf/dzh/dzhi/dzhiq/dzh2i usage matches the solver byte-for-byte. No off-by-one-half anywhere.
• ekh for scalars, ekm for momentum, consistently (:1375, :1381, :1392 vs :1405, :1411, :1421).

Two pre-existing quirks worth noting but not unique to this branch: ekm/ekh include the molecular part, so these outputs are really total diffusive fluxes (and in DNS mode are purely molecular yet still labelled SGS); and the statistics' emom carries dxf/dxhi weighting where the solver uses uniform dxi, identical today but divergent if x-stretching is ever enabled.

9. No test coverage

grep -ri sgs tests/ tools/python/tests/ returns zero matches on this branch — same as master. Nothing exercises any SGS flux, its sign, its staggering, or its IBM masking.

Suggested minimum: a uniform-shear laminar case asserting usgs == -nu_t*du/dz to machine precision would catch items 1, 2 and 6 at once; an IBM case asserting that k=kb of every masked profile is the fill value would catch item 4. The #306 PR adds exactly this kind of in-solver kernel test against master's modstatsdump.f90 — it should be forward-ported to out_stats.f90, ideally by extracting the SGS block into a callable kernel so the test pins the real code rather than a copy of the stencil.