Describe the bug
In idealized simulations, the routine set_physical_bc3d sets the boundary zone points for periodic, symmetric, and open boundary conditions:

At first, the routine does not set the corner points (e.g., (ids-1,jds-1)): For the y-boundaries it iterates i only from MAX(ids,its-1) to MIN(ite+1,ide+istag) and analogously for the x-boundaries (note that for open BC it actually also tries to set the corner points (see line 896), but not effectively, because the RHS values in lines 903-905 are still 0 there, as the y-BC has not been applied yet).

For doubly-periodic BC, the corner points are set correctly afterwards, starting at line 1100. In any other case, however, it seems that the corner points are not assigned and appear as zeros.
At least that's what I see, when debugging the code with gdb.

These corner points are later used, e.g. in the diffusion module:

Thus, because rho(i_start-1,k,j_start-1)=0, xkxavg(i_start, k, j_start) is much smaller than at the other locations.
Inaccuracies must obviously arise.

Suggested solution
The corner points should be set also for non-doubly-periodic boundary conditions. This should be easy to achieve by applying lines 747-1096 in module_bc.F twice, the second time only for the corner points. But someone needs to look at this in more detail.