Vertical integral of vertical viscous tendency doesn't match boundary stress

[clpwolfe]

The vertical integral of the vertical viscous stress should match the stress at the boundaries. I'm running a simple zonal channel in layer mode with bottom drag and no surface stress. I have a flat bottom. The diagnostic `h_du_dt_visc` should be the viscous tendency times the thickness, so the vertical sum of this diagnostic should equal the (negative of the) bottom stress `taux_bot` divided by the Boussinesq density rho0. However, it doesn't. In most places, the relative error between the two is a nearly spatially uniform 2.55e-5, although it gets above 60% at a few isolated grid points. My grid is 1600 x 600, so a few grid points is a vanishingly small part of the domain, but it's weird that there should be such large errors anywhere.

Has anybody else using these diagnostics seen something similar? Anybody know what to make it it? Could it be due to the diagnostic not properly handling layer outcropping?

Thanks!