Looking further, I find that
print *, ubound (soda, 1), ubound (soda, 2)
print *, ubound (soda)
outputs on the first call to the procedure:
0 2
0 0
So, I think something is not right with the vector version of
ubound for this case.
Since I am convinced that this patch is doing the right thing, I
have replaced
if (.not.any(lbound (soda) <= ubound (soda))) call abort ()
with what should be the equivalent
if ((lbound (soda, 1) > ubound (soda, 1)) .and. &
(lbound (soda, 2) > ubound (soda, 2))) call abort ()