bug: Incorrect node provided by "BCobj".stress_period_data.df for some BC packages

[pfjuckem-usgs]

Describe the bug
The "BCobj".stress_period_data.df value for some BC packages returns a 2D node identifier rather than a 3D node identifier. I have only checked the RIV and WEL packages in MF2005 (MF-NWT); other packages and MF6 may or maynot be affected.

To Reproduce
Steps to reproduce the behavior:

1. Unzip the attached zip file.

ex01_mf2005_addriv.zip

2. Run the bad_nodes.py script.
3. The bad_nodes.py script will print (to screen & summary file) dataframes for the WEL and RIV packages from a demo model used by the MODPATH7 documentation. The results will show that the returned "node" column contains values that are accurate only for BC cells in layer 1. For BC cells in deeper layers, the node column values are incorrect. The correct "3D" values are provided in the "node3d" column.
4. Alternatively, download the ex01_mf2005 example from the MODPATH7 repo, add RIV cells to layers deeper layer 1, then run the script on those updated NAM and RIV files. I realize that some may find it odd to add RIV or other surface water BC packages to layers other than layer 1 (I've never done this myself), but I've been "repurposing" several existing models and have encountered numerous instances where SW BC packages include BCs in layers greater than 1, so this definitely happens in practice.

Expected behavior
I expected the "node" values returned by *.stress_period_data.df to match the "node3d" values.

Screenshots
This screen shot shows the output of the summary file.

Desktop (please complete the following information):

• OS: Windows 11
• Browser: N/A
• Version: flopy 3.9.1 (but I didn't see any commits suggesting this had been fixed recently)