Flopy Error When Loading USG Bas6 package with Externally read Starting Heads (STRT)

[javgs-bd]

Hi,

I'm trying to load a mf-usg model that reads initial heads from an external binary head file. The bas6 package has been generated from Groundwater Vistas, and it structured in the followig way:

• The NAM file includes DATA(BINARY) 44 initial_heads.hds stating the unit file for the STRT file will be 44.
• The bas package includes at the end:

EXTERNAL -44  1  (BINARY)  -1  Starting Heads Layer 1
EXTERNAL -44  1  (BINARY)  -1  Starting Heads Layer 2
EXTERNAL -44  1  (BINARY)  -1  Starting Heads Layer 3
EXTERNAL -44  1  (BINARY)  -1  Starting Heads Layer 4
EXTERNAL -44  1  (BINARY)  -1  Starting Heads Layer 5
EXTERNAL -44  1  (BINARY)  -1  Starting Heads Layer 6

Which shows that the STRT array is read from unit 44. Now up to here what puzzles me is that the minus sign preceding the unit file number. I have searched in all the MODFLOW (2000 - MF6) I/O manuals and can not find a specific reference to thye preceding minus sign. The closest idea I got is a ChatGPT answer that states the minus sign is a FORTRAN convention for output files. If I remove that minus sign the reading (U2DREL?) process fails.

The ultimate problem is that the standard flopy flopy.mfusg.MfUsg.load command fails when attempting to load the model, with the message described below. I wonder if there is any source I can be pointed at where I can figure out this problem, and if this is actualy a missing feature in the mfUSG.load method.

I'll appreciate any tips.
Javier

   could not determine filename for unit -44

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
C:\Users\JAVIER~1\AppData\Local\Temp/ipykernel_20336/1887670592.py in <module>
     17 #--- https://stackoverflow.com/questions/51734703/extract-heads-from-modflow-usg-binary-output-using-flopy
     18 
---> 19 gwf = flopy.mfusg.MfUsg.load(
     20                                    model_name,
     21                                    version      = 'mfusg',

~\.conda\envs\flopy-dev\lib\site-packages\flopy\mfusg\mfusg.py in load(cls, f, version, exe_name, verbose, model_ws, load_only, forgive, check)
    261         cls._set_mfpar_pval(model, ext_unit_dict, ext_pkg_d)
    262 
--> 263         files_successfully_loaded, files_not_loaded = cls._load_packages(
    264             model, ext_unit_dict, ext_pkg_d, load_only, forgive
    265         )

~\.conda\envs\flopy-dev\lib\site-packages\flopy\mfusg\mfusg.py in _load_packages(cls, model, ext_unit_dict, ext_pkg_d, load_only, forgive)
    340                     files_successfully_loaded,
    341                     files_not_loaded,
--> 342                 ) = cls._load_ext_unit_dict_paks(
    343                     model,
    344                     ext_unit_dict,

~\.conda\envs\flopy-dev\lib\site-packages\flopy\mfusg\mfusg.py in _load_ext_unit_dict_paks(cls, model, ext_unit_dict, load_only, item, forgive, files_successfully_loaded, files_not_loaded)
    396             if forgive:
    397                 try:
--> 398                     cls._ext_unit_d_load(model, ext_unit_dict, item)
    399                     files_successfully_loaded.append(item.filename)
    400                     if model.verbose:

~\.conda\envs\flopy-dev\lib\site-packages\flopy\mfusg\mfusg.py in _ext_unit_d_load(model, ext_unit_dict, ext_unit_d_item)
    457         package_load_args = getfullargspec(ext_unit_d_item.package.load)[0]
    458         if "check" in package_load_args:
--> 459             ext_unit_d_item.package.load(
    460                 ext_unit_d_item.filehandle,
    461                 model,

~\.conda\envs\flopy-dev\lib\site-packages\flopy\modflow\mfbas.py in load(cls, f, model, ext_unit_dict, check, **kwargs)
    360 
    361         # dataset 4 -- strt
--> 362         strt = Util3d.load(
    363             f, model, (nlay, nrow, ncol), np.float32, "strt", ext_unit_dict
    364         )

~\.conda\envs\flopy-dev\lib\site-packages\flopy\utils\util_array.py in load(cls, f_handle, model, shape, dtype, name, ext_unit_dict, array_format)
    911                 nr = nrow
    912                 nc = ncol
--> 913             u2d = Util2d.load(
    914                 f_handle,
    915                 model,

~\.conda\envs\flopy-dev\lib\site-packages\flopy\utils\util_array.py in load(cls, f_handle, model, shape, dtype, name, ext_unit_dict, array_free_format, array_format)
   2913 
   2914         elif cr_dict["type"] == "external":
-> 2915             ext_unit = ext_unit_dict[cr_dict["nunit"]]
   2916             if ext_unit.filehandle is None:
   2917                 raise OSError(

KeyError: -44

[cnicol-gwlogic]

Hi, not sure this is "the answer" or not, but I find that including the negative sign in the name file works.

[javgs-bd]

It does @cnicol-gwlogic, thanks. Both flopy (loads) and the model keeps running.
I'll ask about that minus sign and report back.

[javgs-bd]

Hi,

I asked Jim Rumbaugh about this, and he tells me that the minus sign tells Modflow it is a binary file. As to why other otuput binary files Nunits in other packages are not negative (such as the hds and cbb files), it is because this minus sign is hardcoded in the modflow bas package reading.

I guess this is something worth catching when reading heads from external files, and wonder if it is a transversal MODFLOW convention or only applies to MF-USG.

Should I create flopy issue from this?

[cnicol-gwlogic]

Hi @javgs-bd I'm pretty confident this would also affect mfnwt, mf2005, etc too. I'm surprised it's not been catered for though. It may be worth checking through the examples / autotests for any cases with external binary input read examples (with negative sign convention), and see what their unit number is in the nam file / how / if they are dealt with. It would be an easy fix to trap for anyway, but maybe check it is not dealt with somehow before raising an issue.