pure subroutine cavity_diff_impl(a, b, active, buf)
!! `a - b` with the cavity missing-value convention.
! assumed-shape-ok: diag fill — fires once per output frame (cadence-bounded).
real(wp), intent(in) :: a(:, :), b(:, :), active(:, :) ! assumed-shape-ok: diag fill — cadence-bounded
real(wp), intent(inout) :: buf(:, :, :) ! assumed-shape-ok: diag fill — cadence-bounded
integer :: i, j, nx, ny
real(wp) :: qnan
nx = min(size(buf, 1), size(a, 1), size(b, 1), size(active, 1))
ny = min(size(buf, 2), size(a, 2), size(b, 2), size(active, 2))
qnan = ieee_value(0.0_wp, ieee_quiet_nan)
do concurrent(j=1:ny, i=1:nx)
if (active(i, j) > 0.5_wp) then
buf(i, j, 1) = a(i, j) - b(i, j)
else
buf(i, j, 1) = qnan
end if
end do
end subroutine cavity_diff_impl