fix(ulstrd): package lists for wel, drn, ghb, etc. were not always re…

…ad correctly

* flopy will now read packages that have lists specified with EXTERNAL
* flopy will now read packages that have lists that use the SFAC keyword
* flopy now has a general list reader flopy.utils.flopy_io.ulstrd patterned after the MODFLOW ULSTRD subroutine
* flopy will now use this ulstrd reader to read lists defined with parameters
* With the m.write_input(), binary external files always had REPLACE tacked onto the end.  This was corrected so that REPLACE is only added when output is set to True.
* Closes modflowpy#683

autotest/t003_test.py

@@ -108,11 +108,11 @@ def test_load_nam_mt_nonexistant_file():
 
 
 if __name__ == '__main__':
-    test_loadoc()
-    test_loadoc_nstpfail()
-    test_load_nam_mf_nonexistant_file()
-    test_load_nam_mt_nonexistant_file()
-    # test_loadoc_lenfail()
-    # test_loadfreyberg()
-    # test_loadoahu()
-    # test_loadtwrip()
+    #test_loadoc()
+    #test_loadoc_nstpfail()
+    #test_load_nam_mf_nonexistant_file()
+    #test_load_nam_mt_nonexistant_file()
+    #test_loadoc_lenfail()
+    #test_loadfreyberg()
+    test_loadoahu()
+    #test_loadtwrip()

autotest/t067_ulstrd.py

@@ -0,0 +1,137 @@
+import os
+import numpy as np
+import flopy
+
+tpth = os.path.join('temp', 't067')
+if not os.path.isdir(tpth):
+    os.makedirs(tpth)
+
+
+def test_ulstrd():
+
+    # Create an original model and then manually modify to use
+    # advanced list reader capabilities
+    ws = tpth
+    nlay = 1
+    nrow = 10
+    ncol = 10
+    nper = 3
+
+    # create the ghbs
+    ghbra = flopy.modflow.ModflowGhb.get_empty(20)
+    l = 0
+    for i in range(nrow):
+        ghbra[l] = (0, i, 0, 1., 100. + i)
+        l += 1
+        ghbra[l] = (0, i, ncol - 1, 1., 200. + i)
+        l += 1
+    ghbspd = {0: ghbra}
+
+    # create the drains
+    drnra = flopy.modflow.ModflowDrn.get_empty(2)
+    drnra[0] = (0, 1, int(ncol / 2), .5, 55.)
+    drnra[1] = (0, 2, int(ncol / 2), .5, 75.)
+    drnspd = {0: drnra}
+
+    # create the wells
+    welra = flopy.modflow.ModflowWel.get_empty(2)
+    welra[0] = (0, 1, 1, -5.)
+    welra[1] = (0, nrow - 3, ncol - 3, -10.)
+    welspd = {0: welra}
+
+    m = flopy.modflow.Modflow(modelname='original', model_ws=ws,
+                              exe_name='mf2005')
+    dis = flopy.modflow.ModflowDis(m, nlay=nlay, nrow=nrow, ncol=ncol,
+                                   nper=nper)
+    bas = flopy.modflow.ModflowBas(m)
+    lpf = flopy.modflow.ModflowLpf(m)
+    ghb = flopy.modflow.ModflowGhb(m, stress_period_data=ghbspd)
+    drn = flopy.modflow.ModflowDrn(m, stress_period_data=drnspd)
+    wel = flopy.modflow.ModflowWel(m, stress_period_data=welspd)
+    pcg = flopy.modflow.ModflowPcg(m)
+    oc = flopy.modflow.ModflowOc(m)
+    m.add_external('original.drn.dat', 71)
+    m.add_external('original.wel.bin', 72, binflag=True, output=False)
+    m.write_input()
+
+    # rewrite ghb
+    fname = os.path.join(ws, 'original.ghb')
+    with open(fname, 'w') as f:
+        f.write('{} {}\n'.format(ghbra.shape[0], 0))
+        for kper in range(nper):
+            f.write('{} {}\n'.format(ghbra.shape[0], 0))
+            f.write('open/close original.ghb.dat\n')
+
+    # write ghb list
+    sfacghb = 5
+    fname = os.path.join(ws, 'original.ghb.dat')
+    with open(fname, 'w') as f:
+        f.write('sfac {}\n'.format(sfacghb))
+        for k, i, j, stage, cond in ghbra:
+            f.write('{} {} {} {} {}\n'.format(k + 1, i + 1, j + 1, stage, cond))
+
+    # rewrite drn
+    fname = os.path.join(ws, 'original.drn')
+    with open(fname, 'w') as f:
+        f.write('{} {}\n'.format(drnra.shape[0], 0))
+        for kper in range(nper):
+            f.write('{} {}\n'.format(drnra.shape[0], 0))
+            f.write('external 71\n')
+
+    # write drn list
+    sfacdrn = 1.5
+    fname = os.path.join(ws, 'original.drn.dat')
+    with open(fname, 'w') as f:
+        for kper in range(nper):
+            f.write('sfac {}\n'.format(sfacdrn))
+            for k, i, j, stage, cond in drnra:
+                f.write(
+                    '{} {} {} {} {}\n'.format(k + 1, i + 1, j + 1, stage, cond))
+
+    # rewrite wel
+    fname = os.path.join(ws, 'original.wel')
+    with open(fname, 'w') as f:
+        f.write('{} {}\n'.format(drnra.shape[0], 0))
+        for kper in range(nper):
+            f.write('{} {}\n'.format(drnra.shape[0], 0))
+            f.write('external 72 (binary)\n')
+
+    # create the wells, but use an all float dtype to write a binary file
+    # use one-based values
+    weldt = np.dtype([('k', '<f4'), ('i', '<f4'), ('j', '<f4'), ('q', '<f4'), ])
+    welra = np.recarray(2, dtype=weldt)
+    welra[0] = (1, 2, 2, -5.)
+    welra[1] = (1, nrow - 2, ncol - 2, -10.)
+    fname = os.path.join(ws, 'original.wel.bin')
+    with open(fname, 'wb') as f:
+        welra.tofile(f)
+        welra.tofile(f)
+        welra.tofile(f)
+
+    # no need to run the model
+    #success, buff = m.run_model(silent=True)
+    #assert success, 'model did not terminate successfully'
+
+    # the m2 model will load all of these external files, possibly using sfac
+    # and just create regular list input files for wel, drn, and ghb
+    fname = 'original.nam'
+    m2 = flopy.modflow.Modflow.load(fname, model_ws=ws, verbose=False)
+    m2.name = 'new'
+    m2.write_input()
+
+    originalghbra = m.ghb.stress_period_data[0].copy()
+    originalghbra['cond'] *= sfacghb
+    assert np.array_equal(originalghbra, m2.ghb.stress_period_data[0])
+
+    originaldrnra = m.drn.stress_period_data[0].copy()
+    originaldrnra['cond'] *= sfacdrn
+    assert np.array_equal(originaldrnra, m2.drn.stress_period_data[0])
+
+    originalwelra = m.wel.stress_period_data[0].copy()
+    assert np.array_equal(originalwelra, m2.wel.stress_period_data[0])
+
+    return
+
+
+if __name__ == '__main__':
+    test_ulstrd()

flopy/modflow/mf.py

@@ -429,14 +429,18 @@ def write_name_file(self):
         f_nam.write('{}'.format(self.get_name_file_entries()))
 
         # write the external files
-        for u, f, b in zip(self.external_units, self.external_fnames,
-                           self.external_binflag):
+        for u, f, b, o in zip(self.external_units, self.external_fnames,
+                              self.external_binflag, self.external_output):
             if u == 0:
                 continue
             # fr = os.path.relpath(f, self.model_ws)
+            replace_text = ''
+            if o:
+                replace_text = 'REPLACE'
             if b:
                 f_nam.write(
-                    'DATA(BINARY)   {0:5d}  '.format(u) + f + ' REPLACE\n')
+                    'DATA(BINARY)   {0:5d}  '.format(u) + f +
+                     replace_text + '\n')
             else:
                 f_nam.write('DATA           {0:5d}  '.format(u) + f + '\n')
 

flopy/modflow/mfchd.py

@@ -192,6 +192,10 @@ def get_default_dtype(structured=True):
                               ("ehead", np.float32)])
         return dtype
 
+    @staticmethod
+    def get_sfac_columns():
+        return ['shead', 'ehead']
+
     @staticmethod
     def load(f, model, nper=None, ext_unit_dict=None):
         """

flopy/modflow/mfdrn.py

@@ -241,6 +241,10 @@ def get_empty(ncells=0, aux_names=None, structured=True, is_drt=False):
             dtype = Package.add_to_dtype(dtype, aux_names, np.float32)
         return create_empty_recarray(ncells, dtype, default_value=-1.0E+10)
 
+    @staticmethod
+    def get_sfac_columns():
+        return ['cond']
+
     @staticmethod
     def load(f, model, nper=None, ext_unit_dict=None, check=True):
         """

flopy/modflow/mfghb.py

@@ -226,6 +226,10 @@ def get_default_dtype(structured=True):
                               ("cond", np.float32)])
         return dtype
 
+    @staticmethod
+    def get_sfac_columns():
+        return ['cond']
+
     @staticmethod
     def load(f, model, nper=None, ext_unit_dict=None, check=True):
         """

flopy/modflow/mfhfb.py

@@ -225,6 +225,10 @@ def get_default_dtype(structured=True):
             assert not structured, 'is there an unstructured HFB???'
         return dtype
 
+    @staticmethod
+    def get_sfac_columns():
+        return ['hydchr']
+
     @staticmethod
     def load(f, model, ext_unit_dict=None):
         """
@@ -328,9 +332,6 @@ def load(f, model, ext_unit_dict=None):
                 iname = 'static'
                 par_dict, current_dict = pak_parms.get(pname)
                 data_dict = current_dict[iname]
-                # print par_dict
-                # print data_dict
-
                 par_current = ModflowHfb.get_empty(par_dict['nlst'])
 
                 #
@@ -344,6 +345,7 @@ def load(f, model, ext_unit_dict=None):
 
                 # fill current parameter data (par_current)
                 for ibnd, t in enumerate(data_dict):
+                    t = tuple(t)
                     par_current[ibnd] = tuple(t[:len(par_current.dtype.names)])
 
                 # convert indices to zero-based

flopy/modflow/mfparbc.py

@@ -5,7 +5,7 @@
 """
 
 import numpy as np
-from ..utils.flopy_io import line_strip
+from ..utils.flopy_io import line_strip, ulstrd
 
 
 class ModflowParBc(object):
@@ -42,7 +42,7 @@ def get(self, fkey):
         return None
 
     @staticmethod
-    def load(f, npar, dt, verbose=False):
+    def load(f, npar, dt, model, ext_unit_dict=None, verbose=False):
         """
         Load bc property parameters from an existing bc package
         that uses list data (e.g. WEL, RIV, etc.).
@@ -101,22 +101,9 @@ def load(f, npar, dt, verbose=False):
                         instnam = t[0].lower()
                     else:
                         instnam = 'static'
-                    bcinst = []
-                    for nw in range(nlst):
-                        line = f.readline()
-                        t = line_strip(line).split()
-                        bnd = []
-                        for jdx in range(nitems):
-                            # if jdx < 3:
-                            if issubclass(dt[jdx].type, np.integer):
-                                # conversion to zero-based occurs in package load method in mbase.
-                                bnd.append(int(t[jdx]))
-                            else:
-                                try:
-                                    bnd.append(float(t[jdx]))
-                                except IndexError:
-                                    bnd.append(1.)
-                        bcinst.append(bnd)
+
+                    ra = np.zeros(nlst, dtype=dt)
+                    bcinst = ulstrd(f, nlst, ra, model, [], ext_unit_dict)
                     pinst[instnam] = bcinst
                 bc_parms[parnam] = [{'partyp': partyp, 'parval': parval,
                                      'nlst': nlst, 'timevarying': timeVarying},

flopy/modflow/mfriv.py

@@ -253,6 +253,10 @@ def get_default_dtype(structured=True):
 
         return dtype
 
+    @staticmethod
+    def get_sfac_columns():
+        return ['cond']
+
     def ncells(self):
         # Return the  maximum number of cells that have river
         # (developed for MT3DMS SSM package)

flopy/modflow/mfwel.py

@@ -331,6 +331,10 @@ def get_empty(ncells=0, aux_names=None, structured=True):
             dtype = Package.add_to_dtype(dtype, aux_names, np.float32)
         return create_empty_recarray(ncells, dtype, default_value=-1.0E+10)
 
+    @staticmethod
+    def get_sfac_columns():
+        return ['flux']
+
     @staticmethod
     def load(f, model, nper=None, ext_unit_dict=None, check=True):
         """