A bounded coordinatewise nondecreasing function on a cube is Riemann integrable

For a positive integer $n$, partition each coordinate into $[k\mathrm{/}n,(k+1)\mathrm{/}n]$, $0\le k<n$, and write $C_k$ for the product cell with lower vertex ${\mathrm{\ell }}_k=(k_1\mathrm{/}n,\dots ,k_d\mathrm{/}n)$ and upper vertex $u_k=((k_1+1)\mathrm{/}n,\dots ,(k_d+1)\mathrm{/}n)$. Coordinatewise monotonicity gives ${\inf }_{C_k}h=h({\mathrm{\ell }}_k)$ and ${\sup }_{C_k}h=h(u_k)$. Each cell volume is $n^{-d}$, so its contribution to $U-L$ is $n^{-d}(h(u_k)-h({\mathrm{\ell }}_k))$.

For each $k$, let $z_k^{(j)}$ have its first $j$ coordinates at the upper knots of $C_k$ and the remaining coordinates at the lower knots, for $0\le j\le d$. Then $z_k^{(0)}={\mathrm{\ell }}_k$, $z_k^{(d)}=u_k$, and $h(u_k)-h({\mathrm{\ell }}_k)={\sum }_{j=1}^d(h(z_k^{(j)})-h(z_k^{(j-1)}))$; all summands are nonnegative. Fix $j$ and the other $d-1$ indices. Summing its increments over $k_j=0,\dots ,n-1$ telescopes to $h(y_1)-h(y_0)$, where $y_0,y_1$ have the other coordinates at the selected lower or upper knots and the $j$th coordinates $0,1$. Since $\mathbf{0}\le y_0\le y_1\le \mathbf{1}$, this is at most $h(\mathbf{1})-h(\mathbf{0})$. There are $n^{d-1}$ choices of the other indices for each $j$, so ${\sum }_k(h(u_k)-h({\mathrm{\ell }}_k))\le d{n}^{d-1}(h(\mathbf{1})-h(\mathbf{0}))$.

By step 1.1, the grid's Darboux gap is at most $d(h(\mathbf{1})-h(\mathbf{0}))\mathrm{/}n$. Choose $n$ large enough using [F4] to make this less than any prescribed positive tolerance (if the range is zero, every gap is zero). The criterion [F3] then gives Riemann integrability.