fix(upw): keyword for vertical K (#746)

autotest/t003_test.py

@@ -44,6 +44,20 @@ def test_loadtwrip():
     assert ml.load_fail is False
 
     return
+
+
+def test_loadtwrip_upw():
+    cwd = os.getcwd()
+    pth = os.path.join('..', 'examples', 'data', 'parameters')
+    assert (os.path.isdir(pth))
+    os.chdir(pth)
+    namefile = 'twrip_upw.nam'
+    ml = flopy.modflow.Modflow.load(namefile, verbose=True)
+    os.chdir(cwd)
+    assert isinstance(ml, flopy.modflow.Modflow)
+    assert ml.load_fail is False
+
+    return
 
 
 def test_loadoc():
@@ -116,3 +130,4 @@ def test_load_nam_mt_nonexistant_file():
     test_loadfreyberg()
     test_loadoahu()
     test_loadtwrip()
+    test_loadtwrip_upw()

examples/data/parameters/twrip_upw.nam

@@ -0,0 +1,14 @@
+LIST   2  twrip_upw.lst
+BAS6   3  twrip.ba6
+PVAL   4  twrip_upw.pval
+UPW   11  twrip_upw.upw
+WEL   12  twrip.wel
+DRN   13  twrip.drn
+RCH   18  twrip.rch
+NWT   19  twrip_upw.nwt
+OC    22  twrip.oc
+MULT   8  twrip.mlt
+ZONE   9  twrip_upw.zon
+DIS   10  twrip.dis
+DATA  101 twrip.ib1.ref
+DATA  102 rchzones.ref

examples/data/parameters/twrip_upw.nwt

@@ -0,0 +1,2 @@
+# NWT package for  MODFLOW-NWT, generated by Flopy.
+ 1.000e-02 5.000e+02       100 1.000e-05         1         0         0    SIMPLE

examples/data/parameters/twrip_upw.pval

@@ -0,0 +1,6 @@
+# Example pval
+4
+HK1_z1    1.0E-3
+HK1_z2    1.0E-2
+VK1_z1    5.0E-4
+VK1_z2    5.0E-3

examples/data/parameters/twrip_upw.upw

@@ -0,0 +1,35 @@
+#  UPW Package input data for example problem
+0  1.00E+30    10    0 IUPWCB HDRY NPUPW IPHDRY
+ 1 0 0
+ 0 0 0
+ 1. 1. 1.
+ 0 0 0
+ 0 0 0
+HK1_z1 HK  1.0E-3 1
+1 NONE upwzones 1
+HK1_z2 HK  1.0E-3 1
+1 NONE upwzones 2
+HK2 HK  1.0E-4 1
+2 NONE ALL
+HK3 HK  2.0E-4 1
+3 NONE ALL
+VK1_z1 VK  1.0E-4 1
+1 NONE upwzones 1
+VK1_z2 VK  1.0E-4 1
+1 NONE upwzones 2
+VK2 VK 1.0E-5 1
+2 NONE ALL
+VK3 VK 1.0E-5 1
+3 NONE ALL
+VKCB1 VKCB 1.0 1
+1 MULT1 ALL
+VKCB2 VKCB 0.5 1
+2 MULT1 ALL
+0                        HK layer   1
+0                        VK layer   1
+0                        VKCB layer 1
+0                        HK layer   2
+0                        VK layer   2
+0                        VKCB layer 2
+0                        HK  ayer   3
+0                        VK layer   3

examples/data/parameters/twrip_upw.zon

@@ -0,0 +1,20 @@
+2
+RCHZONES
+EXTERNAL 102 1 (free) 1
+upwzones
+INTERNAL 1 (15I1) -1
+111111111111111
+111111111111111
+111111111111111
+111111111111111
+111111111111111
+111111111111111
+111111111111111
+111111111111111
+111111111111111
+111111111111111
+111111112211111
+111111122221111
+111111111222111
+111111111122222
+111111111111222

flopy/modflow/mfhob.py

@@ -577,9 +577,9 @@ def __init__(self, model, tomulth=1., obsname='HOBS',
             for key, value in self.mlay.items():
                 tot += value
             if not (np.isclose(tot, 1.0, rtol=0)):
-                msg = 'sum of dataset 4 proportions must equal 1.0 - ' + \
+                msg = ('sum of dataset 4 proportions must equal 1.0 - ' + \
                        'sum of dataset 4 proportions = {tot} for ' + \
-                       'observation name {obsname}.'.format(
+                       'observation name {obsname}.').format(
                            tot=tot,
                            obsname=self.obsname)
                 raise ValueError(msg)

flopy/modflow/mfupw.py

@@ -251,11 +251,12 @@ def write_file(self, check=True, f=None):
         # Item 0: text
         f_upw.write('{}\n'.format(self.heading))
         # Item 1: IBCFCB, HDRY, NPLPF
-        f_upw.write('{0:10d}{1:10.3G}{2:10d}{3:10d}{4:s}\n'.format(self.ipakcb,
-                                                                   self.hdry,
-                                                                   self.npupw,
-                                                                   self.iphdry,
-                                                                   self.options))
+        f_upw.write('{0:10d}{1:10.3G}{2:10d}{3:10d}{4:s}\n'
+                    .format(self.ipakcb,
+                            self.hdry,
+                            self.npupw,
+                            self.iphdry,
+                            self.options))
         # LAYTYP array
         f_upw.write(self.laytyp.string)
         # LAYAVG array
@@ -276,7 +277,7 @@ def write_file(self, check=True, f=None):
             if self.chani[k] < 1:
                 f_upw.write(self.hani[k].get_file_entry())
             f_upw.write(self.vka[k].get_file_entry())
-            if transient == True:
+            if transient:
                 f_upw.write(self.ss[k].get_file_entry())
                 if self.laytyp[k] != 0:
                     f_upw.write(self.sy[k].get_file_entry())
@@ -389,8 +390,7 @@ def load(f, model, ext_unit_dict=None, check=True):
         laywet = np.empty((nlay,), dtype=np.int32)
         laywet = read1d(f, laywet)
 
-        # Item 7: WETFCT, IWETIT, IHDWET
-        wetfct, iwetit, ihdwet = None, None, None
+        # check that LAYWET is 0 for all layers
         iwetdry = laywet.sum()
         if iwetdry > 0:
             raise Exception('LAYWET should be 0 for UPW')
@@ -408,7 +408,10 @@ def load(f, model, ext_unit_dict=None, check=True):
         ss = [0] * nlay
         sy = [0] * nlay
         vkcb = [0] * nlay
+        # load by layer
         for k in range(nlay):
+
+            # hk
             if model.verbose:
                 print('   loading hk layer {0:3d}...'.format(k + 1))
             if 'hk' not in par_types:
@@ -419,6 +422,8 @@ def load(f, model, ext_unit_dict=None, check=True):
                 t = mfpar.parameter_fill(model, (nrow, ncol), 'hk', parm_dict,
                                          findlayer=k)
             hk[k] = t
+
+            # hani
             if chani[k] < 1:
                 if model.verbose:
                     print('   loading hani layer {0:3d}...'.format(k + 1))
@@ -430,23 +435,26 @@ def load(f, model, ext_unit_dict=None, check=True):
                     t = mfpar.parameter_fill(model, (nrow, ncol), 'hani',
                                              parm_dict, findlayer=k)
                 hani[k] = t
+
+            # vka
             if model.verbose:
                 print('   loading vka layer {0:3d}...'.format(k + 1))
+            key = 'vk'
+            if layvka[k] != 0:
+                key = 'vani'
             if 'vk' not in par_types and 'vani' not in par_types:
-                key = 'vka'
-                if layvka[k] != 0:
-                    key = 'vani'
                 t = Util2d.load(f, model, (nrow, ncol), np.float32, key,
                                 ext_unit_dict)
             else:
                 line = f.readline()
-                key = 'vka'
-                if 'vani' in par_types:
-                    key = 'vani'
                 t = mfpar.parameter_fill(model, (nrow, ncol), key, parm_dict,
                                          findlayer=k)
             vka[k] = t
+
+            # storage properties
             if transient:
+
+                # ss
                 if model.verbose:
                     print('   loading ss layer {0:3d}...'.format(k + 1))
                 if 'ss' not in par_types:
@@ -457,6 +465,8 @@ def load(f, model, ext_unit_dict=None, check=True):
                     t = mfpar.parameter_fill(model, (nrow, ncol), 'ss',
                                              parm_dict, findlayer=k)
                 ss[k] = t
+
+                # sy
                 if laytyp[k] != 0:
                     if model.verbose:
                         print('   loading sy layer {0:3d}...'.format(k + 1))
@@ -469,6 +479,8 @@ def load(f, model, ext_unit_dict=None, check=True):
                         t = mfpar.parameter_fill(model, (nrow, ncol), 'sy',
                                                  parm_dict, findlayer=k)
                     sy[k] = t
+
+            # vkcb
             if model.get_package('DIS').laycbd[k] > 0:
                 if model.verbose:
                     print('   loading vkcb layer {0:3d}...'.format(k + 1))