Mask topsystem packages where IBOUND == -1

In issue #1766 it turned out that iMOD5 filters away RIV & DRN cells where IBOUND == -1 (and thus CHD cells are created), as these create shortcircuits in model cells, potentially leading to unnecessary stability issues.

We can approach this in two ways:

• Only implement this in the from_imod5 logic. We can then relatively easily use IBOUND == -1 as mask.
• Generalize this: Create a method on the groundwatermodel flow object to filter all GHB, RIV, & DRN cells where a CHD cell is located. MODFLOW6 only accepts one CHD package, so we don't require any logic to merge CHD packages together for this masking. This can then also be called when calling Modflow6Simulation.clip_box with states_for_boundary.

I prefer the latter option.

Start with the from_imod5 logic and the groundwatermodel flow object, then inspect Modflow6Simulation.clip_box with states_for_boundary. Determine how GHB, RIV, and DRN cells are filtered when a CHD cell occupies the same location, and consider the preferred generalized behavior. Done means IBOUND == -1 and overlapping CHD cells no longer leave conflicting boundary cells in either entry point.