"""
mfsub module.  Contains the ModflowSub class. Note that the user can access
the ModflowSub class as `flopy.modflow.ModflowSub`.

Additional information for this MODFLOW package can be found at the `Online
MODFLOW Guide
<http://water.usgs.gov/ogw/modflow/MODFLOW-2005-Guide/sub.htm>`_.

"""
import sys

import numpy as np

from ..pakbase import Package
from ..utils import Util2d, Util3d, read1d


class ModflowSub(Package):
    """
    MODFLOW SUB Package Class.

    Parameters
    ----------
    model : model object
        The model object (of type :class:`flopy.modflow.mf.Modflow`) to which
        this package will be added.
    ipakcb : int
        A flag that is used to determine if cell-by-cell budget data should be
        saved. If ipakcb is non-zero cell-by-cell budget data will be saved.
        (default is 0).
    isuboc : int
        isuboc is a flag used to control output of information generated by the
        SUB Package. (default is 0).
    idsave : int
        idsave is a flag and a unit number on which restart records for delay
        interbeds will be saved at the end of the simulation. (default is 0).
    idrest : int
        idrest is a flag and a unit number on which restart records for delay
        interbeds will be read in at the start of the simulation (default is 0).
    idbit : int
        idrest is an optional flag that defines if iteration will be used for
        the delay bed solution of heads. If idbit is 0 or less than 0,
        iteration will not be used (this is is identical to the approach used
        in MODFLOW-2005 versions less than 1.13. if idbit is greater than
        0, iteration of the delay bed solution will continue until convergence
        is achieved. (default is 0).
    nndb : int
        nndb is the number of systems of no-delay interbeds. (default is 1).
    ndb : int
        ndb is the number of systems of delay interbeds. (default is 1).
    nmz : int
        nmz is the number of material zones that are needed to define the
        hydraulic properties of systems of delay interbeds. Each material zone
        is defined by a combination of vertical hydraulic conductivity, elastic
        specific storage, and inelastic specific storage. (default is 1).
    nn : int
        nn is the number of nodes used to discretize the half space to
        approximate the head distributions in systems of delay interbeds.
        (default is 20).
    ac1 : float
        ac1 is an acceleration parameter.  This parameter is used to predict
        the aquifer head at the interbed boundaries on the basis of the head
        change computed for the previous iteration. A value of 0.0 results in
        the use of the aquifer head at the previous iteration. Limited
        experience indicates that optimum values may range from 0.0 to 0.6.
        (default is 0).
    ac2 : float
        ac2 is an acceleration parameter. This acceleration parameter is a
        multiplier for the head changes to compute the head at the new
        iteration. Values are normally between 1.0 and 2.0, but the optimum is
        probably closer to 1.0 than to 2.0. However this parameter also can be
        used to help convergence of the iterative solution by using values
        between 0 and 1. (default is 1.0).
    itmin : int
        ITMIN is the minimum number of iterations for which one-dimensional
        equations will be solved for flow in interbeds when the Strongly
        Implicit Procedure (SIP) is used to solve the ground-water flow
        equations. If the current iteration level is greater than ITMIN and
        the SIP convergence criterion for head closure (HCLOSE) is met at a
        particular cell, the one-dimensional equations for that cell will not
        be solved. The previous solution will be used. The value of ITMIN is
        not used if a solver other than SIP is used to solve the ground-water
        flow equations. (default is 5).
    ln : int or array of ints (nndb)
        ln is a one-dimensional array specifying the model layer assignments
        for each system of no-delay interbeds. (default is 0).
    ldn : int or array of ints (ndb)
        ldn is a one-dimensional array specifying the model layer assignments
        for each system of delay interbeds.(default is 0).
    rnb : float or array of floats (ndb, nrow, ncol)
        rnb is an array specifying the factor nequiv at each cell for each
        system of delay interbeds. The array also is used to define the areal
        extent of each system of interbeds. For cells beyond the areal extent
        of the system of interbeds, enter a number less than 1.0 in the
        corresponding element of this array. (default is 1).
    hc : float or array of floats (nndb, nrow, ncol)
        hc is an array specifying the preconsolidation head or preconsolidation
        stress in terms of head in the aquifer for systems of no-delay
        interbeds. For any model cells in which specified HC is greater than
        the corresponding value of starting head, the value of HC will be set
        to that of starting head. (default is 100000).
    sfe : float or array of floats (nndb, nrow, ncol)
        sfe is an array specifying the dimensionless elastic skeletal storage
        coefficient for systems of no-delay interbeds. (default is 1.e-4).
    sfv : float or array of floats (nndb, nrow, ncol)
        sfv is an array specifying the dimensionless inelastic skeletal storage
        coefficient for systems of no-delay interbeds. (default is 1.e-3).
    com : float or array of floats (nndb, nrow, ncol)
        com is an array specifying the starting compaction in each system of
        no-delay interbeds. Compaction values computed by the package are added
        to values in this array so that printed or stored values of compaction
        and land subsidence may include previous components. Values in this
        array do not affect calculations of storage changes or resulting
        compaction. For simulations in which output values are to reflect
        compaction and subsidence since the start of the simulation, enter zero
        values for all elements of this array. (default is 0).
    dp : list or array of floats (nmz, 3)
        Data item includes nmz records, each with a value of vertical hydraulic
        conductivity, elastic specific storage, and inelastic specific storage.
        (default is [1.e-6, 6.e-6, 6.e-4]).
    dstart : float or array of floats (ndb, nrow, ncol)
        dstart is an array specifying starting head in interbeds for systems of
        delay interbeds. For a particular location in a system of interbeds,
        the starting head is applied to every node in the string of nodes that
        approximates flow in half of a doubly draining interbed.
        (default is 1).
    dhc : float or array of floats (ndb, nrow, ncol)
        dhc is an array specifying the starting preconsolidation head in
        interbeds for systems of delay interbeds. For a particular location in
        a system of interbeds, the starting preconsolidation head is applied to
        every node in the string of nodes that approximates flow in half of a
        doubly draining interbed. For any location at which specified starting
        preconsolidation head is greater than the corresponding value of the
        starting head, Dstart, the value of the starting preconsolidation head
        will be set to that of the starting head. (default is 100000).
    dcom : float or array of floats (ndb, nrow, ncol)
        dcom is an array specifying the starting compaction in each system of
        delay interbeds. Compaction values computed by the package are added to
        values in this array so that printed or stored values of compaction and
        land subsidence may include previous components. Values in this array
        do not affect calculations of storage changes or resulting compaction.
        For simulations in which output values are to reflect compaction
        and subsidence since the start of the simulation, enter zero values
        for all elements of this array. (default is 0).
    dz : float or array of floats (ndb, nrow, ncol)
        dz is an array specifying the equivalent thickness for a system of
        delay interbeds. (default is 1).
    nz : int or array of ints (ndb, nrow, ncol)
        nz is an array specifying the material zone numbers for systems of
        delay interbeds. The zone number for each location in the model grid
        selects the hydraulic conductivity, elastic specific storage, and
        inelastic specific storage of the interbeds. (default is 1).
    ids15 : list or array of ints (12)
        Format codes and unit numbers for subsidence, compaction by model
        layer, compaction by interbed system, vertical displacement, no-delay
        preconsolidation, and delay preconsolidation will be printed. If ids15
        is None and isuboc>0 then print code 0 will be used for all data which
        is output to the binary subsidence output file (unit=1051). The 12
        entries in ids15 correspond to ifm1, iun1, ifm2, iun2, ifm3, iun3,
        ifm4, iun4, ifm5, iun5, ifm6, and iun6 variables. (default is None).
    ids16 : list or array of ints (isuboc, 17)
        Stress period and time step range and print and save flags used to
        control printing and saving of information generated by the SUB Package
        during program execution. Each row of ids16 corresponds to isp1, isp2,
        its1, its2, ifl1, ifl2, ifl3, ifl4, ifl5, ifl6, ifl7, ifl8, ifl9,
        ifl10, ifl11, ifl12, and ifl13 variables for isuboc entries. isp1,
        isp2, its1, and its2 are stress period and time step ranges. ifl1 and
        ifl2 control subsidence printing and saving. ifl3 and ifl4 control
        compaction by model layer printing and saving. ifl5 and ifl6 control
        compaction by interbed system printing and saving. ifl7 and ifl8
        control vertical displacement printing and saving. ifl9 and ifl10
        control critical head for no-delay interbeds printing and saving.
        ifl11 and ifl12 control critical head for delay interbeds printing
        and saving. ifl13 controls volumetric budget for delay interbeds
        printing. If ids16 is None and isuboc>0 then all available subsidence
        output will be printed and saved to the binary subsidence output file
        (unit=1051). (default is None).
    unitnumber : int
        File unit number (default is None).
    filenames : str or list of str
        Filenames to use for the package and the output files. If
        filenames=None the package name will be created using the model name
        and package extension and the cbc output name and other sub output
        files will be created using the model name and .cbc and swt output
        extensions (for example, modflowtest.cbc), if ipakcbc and other
        sub output files (dataset 15) are numbers greater than zero.
        If a single string is passed the package name will be set to the
        string and other sub output files will be set to the model name with
        the appropriate output file extensions. To define the names for all
        package files (input and output) the length of the list of strings
        should be 9. Default is None.

    Attributes
    ----------

    Methods
    -------

    See Also
    --------

    Notes
    -----
    Parameters are supported in Flopy only when reading in existing models.
    Parameter values are converted to native values in Flopy and the
    connection to "parameters" is thus nonexistent. Parameters are not
    supported in the SUB Package.

    Examples
    --------

    >>> import flopy
    >>> m = flopy.modflow.Modflow()
    >>> sub = flopy.modflow.ModflowSub(m)

    """

    def __init__(self, model, ipakcb=None, isuboc=0, idsave=None, idrest=None,
                 idbit=None,
                 nndb=1, ndb=1, nmz=1, nn=20, ac1=0., ac2=1.0, itmin=5,
                 ln=0, ldn=0, rnb=1,
                 hc=100000., sfe=1.e-4, sfv=1.e-3, com=0.,
                 dp=[[1.e-6, 6.e-6, 6.e-4]],
                 dstart=1., dhc=100000., dcom=0., dz=1., nz=1,
                 ids15=None, ids16=None,
                 extension='sub', unitnumber=None,
                 filenames=None):
        """
        Package constructor.

        """
        # set default unit number of one is not specified
        if unitnumber is None:
            unitnumber = ModflowSub.defaultunit()

        # set filenames
        if filenames is None:
            filenames = [None for x in range(9)]
        elif isinstance(filenames, str):
            filenames = [filenames] + [None for x in range(8)]
        elif isinstance(filenames, list):
            if len(filenames) < 9:
                n = 9 - len(filenames) + 1
                filenames = filenames + [None for x in range(n)]

        # update external file information with cbc output, if necessary
        if ipakcb is not None:
            fname = filenames[1]
            model.add_output_file(ipakcb, fname=fname,
                                  package=ModflowSub.ftype())
        else:
            ipakcb = 0

        if idsave is not None:
            fname = filenames[2]
            model.add_output_file(idsave, fname=fname, extension='rst',
                                  package=ModflowSub.ftype())
        else:
            idsave = 0

        if idrest is None:
            idrest = 0

        item15_extensions = ["subsidence.hds", "total_comp.hds",
                             "inter_comp.hds", "vert_disp.hds",
                             "nodelay_precon.hds", "delay_precon.hds"]
        item15_units = [2052 + i for i in range(len(item15_extensions))]

        if isuboc > 0:
            idx = 0
            for k in range(1, 12, 2):
                ext = item15_extensions[idx]
                if ids15 is None:
                    iu = item15_units[idx]
                else:
                    iu = ids15[k]
                fname = filenames[idx+3]
                model.add_output_file(iu, fname=fname, extension=ext,
                                      package=ModflowSub.ftype())
                idx += 1

        extensions = [extension]
        name = [ModflowSub.ftype()]
        units = [unitnumber]
        extra = ['']

        # set package name
        fname = [filenames[0]]

        # Call ancestor's init to set self.parent, extension, name and
        # unit number
        Package.__init__(self, model, extension=extensions, name=name,
                         unit_number=units, extra=extra, filenames=fname)

        nrow, ncol, nlay, nper = self.parent.nrow_ncol_nlay_nper

        self.heading = '# {} package for '.format(self.name[0]) + \
                       ' {}, '.format(model.version_types[model.version]) + \
                       'generated by Flopy.'
        self.url = 'sub.htm'

        self.ipakcb = ipakcb
        self.isuboc = isuboc
        self.idsave = idsave
        self.idrest = idrest
        self.idbit = idbit
        self.nndb = nndb
        self.ndb = ndb
        self.nmz = nmz
        self.nn = nn
        self.ac1 = ac1
        self.ac2 = ac2
        self.itmin = itmin
        # no-delay bed data
        self.ln = None
        self.hc = None
        self.sfe = None
        self.sfv = None
        if nndb > 0:
            self.ln = Util2d(model, (nndb,), np.int32, ln, name='ln')
            self.hc = Util3d(model, (nndb, nrow, ncol), np.float32, hc,
                             name='hc',
                             locat=self.unit_number[0])
            self.sfe = Util3d(model, (nndb, nrow, ncol), np.float32, sfe,
                              name='sfe',
                              locat=self.unit_number[0])
            self.sfv = Util3d(model, (nndb, nrow, ncol), np.float32, sfv,
                              name='sfv',
                              locat=self.unit_number[0])
            self.com = Util3d(model, (nndb, nrow, ncol), np.float32, com,
                              name='com',
                              locat=self.unit_number[0])
        # delay bed data
        self.ldn = None
        self.rnb = None
        self.dstart = None
        self.dhc = None
        self.dz = None
        self.nz = None
        if ndb > 0:
            self.ldn = Util2d(model, (ndb,), np.int32, ldn, name='ldn')
            self.rnb = Util3d(model, (ndb, nrow, ncol), np.float32, rnb,
                              name='rnb',
                              locat=self.unit_number[0])
            self.dstart = Util3d(model, (ndb, nrow, ncol), np.float32, dstart,
                                 name='dstart',
                                 locat=self.unit_number[0])
            self.dhc = Util3d(model, (ndb, nrow, ncol), np.float32, dhc,
                              name='dhc',
                              locat=self.unit_number[0])
            self.dcom = Util3d(model, (ndb, nrow, ncol), np.float32, dcom,
                               name='dcom',
                               locat=self.unit_number[0])
            self.dz = Util3d(model, (ndb, nrow, ncol), np.float32, dz,
                             name='dz',
                             locat=self.unit_number[0])
            self.nz = Util3d(model, (ndb, nrow, ncol), np.int32, nz, name='nz',
                             locat=self.unit_number[0])
        # material zone data
        if isinstance(dp, list):
            dp = np.array(dp)
        self.dp = dp

        # output data
        if isuboc > 0:
            if ids15 is None:
                ids15 = np.zeros(12, dtype=np.int32)
                iu = 0
                for i in range(1, 12, 2):
                    ids15[i] = item15_units[iu]
                    iu += 1
                self.ids15 = ids15
            else:
                if isinstance(ids15, list):
                    ids15 = np.array(ids15)
                self.ids15 = ids15

            if ids16 is None:
                self.isuboc = 1
                # save and print everything
                ids16 = np.ones((1, 17), dtype=np.int32)
                ids16[0, 0] = 0
                ids16[0, 1] = nper - 1
                ids16[0, 2] = 0
                ids16[0, 3] = 9999
            else:
                if isinstance(ids16, list):
                    ids16 = np.array(ids16)
            if len(ids16.shape) == 1:
                ids16 = np.reshape(ids16, (1, ids16.shape[0]))
            self.ids16 = ids16

        # add package to model
        self.parent.add_package(self)

    def write_file(self, check=False, f=None):
        """
        Write the package file.

        Returns
        -------
        None

        """
        if check:
            print("warning: check not implemented for sub")
        nrow, ncol, nlay, nper = self.parent.nrow_ncol_nlay_nper
        # Open file for writing
        if f is None:
            f = open(self.fn_path, 'w')
        # First line: heading
        f.write('{}\n'.format(self.heading))
        # write dataset 1
        f.write(
            '{} {} {} {} {} {} '.format(self.ipakcb, self.isuboc, self.nndb,
                                        self.ndb, self.nmz, self.nn))

        f.write('{} {} {} {} {}'.format(self.ac1, self.ac2,
                                        self.itmin, self.idsave,
                                        self.idrest))
        line = ''
        if self.idbit is not None:
            line += ' {}'.format(self.idbit)
        line += '\n'
        f.write(line)

        if self.nndb > 0:
            t = self.ln.array
            for tt in t:
                f.write('{} '.format(tt + 1))
            f.write('\n')
        if self.ndb > 0:
            t = self.ldn.array
            for tt in t:
                f.write('{} '.format(tt + 1))
            f.write('\n')

        # write dataset 4
        if self.ndb > 0:
            for k in range(self.ndb):
                f.write(self.rnb[k].get_file_entry())

        # write dataset 5 to 8
        if self.nndb > 0:
            for k in range(self.nndb):
                f.write(self.hc[k].get_file_entry())
                f.write(self.sfe[k].get_file_entry())
                f.write(self.sfv[k].get_file_entry())
                f.write(self.com[k].get_file_entry())

        # write dataset 9
        if self.ndb > 0:
            for k in range(self.nmz):
                line = '{:15.6g} {:15.6g} {:15.6g}'.format(self.dp[k, 0],
                                                           self.dp[k, 1],
                                                           self.dp[k, 2]) + \
                       '    #material zone {} data\n'.format(k + 1)
                f.write(line)
        # write dataset 10 to 14
        if self.ndb > 0:
            for k in range(self.ndb):
                f.write(self.dstart[k].get_file_entry())
                f.write(self.dhc[k].get_file_entry())
                f.write(self.dcom[k].get_file_entry())
                f.write(self.dz[k].get_file_entry())
                f.write(self.nz[k].get_file_entry())

        # write dataset 15 and 16
        if self.isuboc > 0:
            # dataset 15
            for i in self.ids15:
                f.write('{} '.format(i))
            f.write('  #dataset 15\n')

            # dataset 16
            for k in range(self.isuboc):
                t = self.ids16[k, :]
                t[0:4] += 1
                for i in t:
                    f.write('{} '.format(i))
                f.write('  #dataset 16 isuboc {}\n'.format(k + 1))

        # close sub file
        f.close()

    @staticmethod
    def load(f, model, ext_unit_dict=None):
        """
        Load an existing package.

        Parameters
        ----------
        f : filename or file handle
            File to load.
        model : model object
            The model object (of type :class:`flopy.modflow.mf.Modflow`) to
            which this package will be added.
        ext_unit_dict : dictionary, optional
            If the arrays in the file are specified using EXTERNAL,
            or older style array control records, then `f` should be a file
            handle.  In this case ext_unit_dict is required, which can be
            constructed using the function
            :class:`flopy.utils.mfreadnam.parsenamefile`.

        Returns
        -------
        sub : ModflowSub object

        Examples
        --------

        >>> import flopy
        >>> m = flopy.modflow.Modflow()
        >>> sub = flopy.modflow.ModflowSub.load('test.sub', m)

        """

        if model.verbose:
            sys.stdout.write('loading sub package file...\n')

        openfile = not hasattr(f, 'read')
        if openfile:
            filename = f
            f = open(filename, 'r')

        # dataset 0 -- header
        while True:
            line = f.readline()
            if line[0] != '#':
                break
        # determine problem dimensions
        nrow, ncol, nlay, nper = model.get_nrow_ncol_nlay_nper()

        # read dataset 1
        if model.verbose:
            sys.stdout.write('  loading sub dataset 1\n')
        t = line.strip().split()
        ipakcb, isuboc, nndb, ndb, nmz, nn = int(t[0]), int(t[1]), int(t[2]), \
                                             int(t[3]), int(t[4]), int(t[5])
        ac1, ac2 = float(t[6]), float(t[7])
        itmin, idsave, idrest = int(t[8]), int(t[9]), int(t[10])

        idbit = None
        if len(t) > 11:
            if isinstance(t[11], (int, float)):
                idbit = int(t[11])
        if idbit is None:
            if model.verbose:
                print('   explicit idbit in file')

        ln = None
        if nndb > 0:
            if model.verbose:
                sys.stdout.write('  loading sub dataset 2\n')
            ln = np.empty((nndb), dtype=np.int32)
            ln = read1d(f, ln) - 1
        ldn = None
        if ndb > 0:
            if model.verbose:
                sys.stdout.write('  loading sub dataset 3\n')
            ldn = np.empty((ndb), dtype=np.int32)
            ldn = read1d(f, ldn) - 1
        rnb = None
        if ndb > 0:
            if model.verbose:
                sys.stdout.write('  loading sub dataset 4\n')
            rnb = [0] * ndb
            for k in range(ndb):
                t = Util2d.load(f, model, (nrow, ncol), np.float32,
                                'rnb delay bed {}'.format(k + 1),
                                ext_unit_dict)
                rnb[k] = t
        hc = None
        sfe = None
        sfv = None
        com = None
        if nndb > 0:
            hc = [0] * nndb
            sfe = [0] * nndb
            sfv = [0] * nndb
            com = [0] * nndb
            for k in range(nndb):
                kk = ln[k] + 1
                # hc
                if model.verbose:
                    sys.stdout.write(
                        '  loading sub dataset 5 for layer {}\n'.format(kk))
                t = Util2d.load(f, model, (nrow, ncol), np.float32,
                                'hc layer {}'.format(kk), ext_unit_dict)
                hc[k] = t
                # sfe
                if model.verbose:
                    sys.stdout.write(
                        '  loading sub dataset 6 for layer {}\n'.format(kk))
                t = Util2d.load(f, model, (nrow, ncol), np.float32,
                                'sfe layer {}'.format(kk), ext_unit_dict)
                sfe[k] = t
                # sfv
                if model.verbose:
                    sys.stdout.write(
                        '  loading sub dataset 7 for layer {}\n'.format(kk))
                t = Util2d.load(f, model, (nrow, ncol), np.float32,
                                'sfv layer {}'.format(kk), ext_unit_dict)
                sfv[k] = t
                # com
                if model.verbose:
                    sys.stdout.write(
                        '  loading sub dataset 8 for layer {}\n'.format(kk))
                t = Util2d.load(f, model, (nrow, ncol), np.float32,
                                'com layer {}'.format(kk), ext_unit_dict)
                com[k] = t

        # dp
        dp = None
        if ndb > 0:
            dp = np.zeros((nmz, 3), dtype=np.float32)
            for k in range(nmz):
                if model.verbose:
                    msg = 2 * ' ' + 'loading sub dataset 9 for material ' + \
                          'zone {}\n'.format(k + 1)
                    sys.stdout.write(msg)
                line = f.readline()
                t = line.strip().split()
                dp[k, :] = float(t[0]), float(t[1]), float(t[2])

        dstart = None
        dhc = None
        dcom = None
        dz = None
        nz = None
        if ndb > 0:
            dstart = [0] * ndb
            dhc = [0] * ndb
            dcom = [0] * ndb
            dz = [0] * ndb
            nz = [0] * ndb
            for k in range(ndb):
                kk = ldn[k] + 1
                # dstart
                if model.verbose:
                    msg = 2 * ' ' + 'loading sub dataset 10 for ' + \
                          'layer {}\n'.format(kk)
                    sys.stdout.write(msg)
                t = Util2d.load(f, model, (nrow, ncol), np.float32,
                                'dstart layer {}'.format(kk),
                                ext_unit_dict)
                dstart[k] = t
                # dhc
                if model.verbose:
                    msg = 2 * ' ' + 'loading sub dataset 11 for ' + \
                          'layer {}\n'.format(kk)
                    sys.stdout.write(msg)
                t = Util2d.load(f, model, (nrow, ncol), np.float32,
                                'dhc layer {}'.format(kk),
                                ext_unit_dict)
                dhc[k] = t
                # dcom
                if model.verbose:
                    msg = 2 * ' ' + 'loading sub dataset 12 for ' + \
                          'layer {}\n'.format(kk)
                    sys.stdout.write(msg)
                t = Util2d.load(f, model, (nrow, ncol), np.float32,
                                'dcom layer {}'.format(kk),
                                ext_unit_dict)
                dcom[k] = t
                # dz
                if model.verbose:
                    msg = 2 * ' ' + 'loading sub dataset 13 for ' + \
                          'layer {}\n'.format(kk)
                    sys.stdout.write(msg)
                t = Util2d.load(f, model, (nrow, ncol), np.float32,
                                'dz layer {}'.format(kk),
                                ext_unit_dict)
                dz[k] = t
                # nz
                if model.verbose:
                    msg = 2 * ' ' + 'loading sub dataset 14 for ' + \
                          'layer {}\n'.format(kk)
                    sys.stdout.write(msg)
                t = Util2d.load(f, model, (nrow, ncol), np.int32,
                                'nz layer {}'.format(kk),
                                ext_unit_dict)
                nz[k] = t

        ids15 = None
        ids16 = None
        if isuboc > 0:
            # dataset 15
            if model.verbose:
                msg = 2 * ' ' + 'loading sub dataset 15 for ' + \
                      'layer {}\n'.format(kk)
                sys.stdout.write(msg)
            ids15 = np.empty(12, dtype=np.int32)
            ids15 = read1d(f, ids15)
            # dataset 16
            ids16 = [0] * isuboc
            for k in range(isuboc):
                if model.verbose:
                    msg = 2 * ' ' + 'loading sub dataset 16 for ' + \
                          'isuboc {}\n'.format(k + 1)
                    sys.stdout.write(msg)
                t = np.empty(17, dtype=np.int32)
                t = read1d(f, t)
                t[0:4] -= 1
                ids16[k] = t

        if openfile:
            f.close()

        model.add_pop_key_list(2051)

        # determine specified unit number
        unitnumber = None
        filenames = [None for x in range(9)]
        if ext_unit_dict is not None:
            unitnumber, filenames[0] = \
                model.get_ext_dict_attr(ext_unit_dict,
                                        filetype=ModflowSub.ftype())
            if ipakcb > 0:
                iu, filenames[1] = \
                    model.get_ext_dict_attr(ext_unit_dict, unit=ipakcb)

            if idsave > 0:
                iu, filenames[2] = \
                    model.get_ext_dict_attr(ext_unit_dict, unit=idsave)

            if isuboc > 0:
                ipos = 3
                for k in range(1, 12, 2):
                    unit = ids15[k]
                    if unit > 0:
                        iu, filenames[ipos] = \
                            model.get_ext_dict_attr(ext_unit_dict, unit=unit)
                        model.add_pop_key_list(unit)
                    ipos += 1

        # create sub instance
        sub = ModflowSub(model, ipakcb=ipakcb, isuboc=isuboc, idsave=idsave,
                         idrest=idrest, idbit=idbit,
                         nndb=nndb, ndb=ndb, nmz=nmz, nn=nn, ac1=ac1, ac2=ac2,
                         itmin=itmin,
                         ln=ln, ldn=ldn, rnb=rnb,
                         hc=hc, sfe=sfe, sfv=sfv, com=com, dp=dp,
                         dstart=dstart, dhc=dhc, dcom=dcom, dz=dz, nz=nz,
                         ids15=ids15, ids16=ids16, unitnumber=unitnumber,
                         filenames=filenames)
        # return sub instance
        return sub

    @staticmethod
    def ftype():
        return 'SUB'

    @staticmethod
    def defaultunit():
        return 32
