How to use the flopy.utils.Util3d function in flopy

the name of the attribute for MfList type.  Only returned if
            attr is MfList.  For example, if attr is MfList and pak is
            flopy.modflow.ModflowWel, then attrname can only be "flux"

        """

        raw = pakattr.lower().split('.')
        if len(raw) != 2:
            self.logger.lraise("pakattr is wrong:{0}".format(pakattr))
        pakname = raw[0]
        attrname = raw[1]
        pak = self.m.get_package(pakname)
        if pak is None:
            if pakname == "extra":
                self.logger.statement("'extra' pak detected:{0}".format(pakattr))
                ud = flopy.utils.Util3d(self.m,(self.m.nlay,self.m.nrow,self.m.ncol),np.float32,1.0,attrname)
                return "extra",ud

            self.logger.lraise("pak {0} not found".format(pakname))
        if hasattr(pak,attrname):
            attr = getattr(pak,attrname)
            return pak,attr
        elif hasattr(pak,"stress_period_data"):
            dtype = pak.stress_period_data.dtype
            if attrname not in dtype.names:
                self.logger.lraise("attr {0} not found in dtype.names for {1}.stress_period_data".\
                                  format(attrname,pakname))
            attr = pak.stress_period_data
            return pak,attr,attrname
        else:
            self.logger.lraise("unrecognized attr:{0}".format(attrname))

self.ict = ict
        self.cinact = cinact
        self.ciclose = ciclose
        self.idisp = idisp
        self.ixdisp = ixdisp
        self.diffnc = diffnc
        self.izod = izod
        self.ifod = ifod
        self.icbund = Util3d(model, (nlay, nrow, ncol), np.float32, icbund,
                             'icbund', )
        self.porosity = Util3d(model, (nlay, nrow, ncol), np.float32,
                               porosity, 'porosity')
        # self.arad = Util2d(model, (1, nja), np.float32,
        #                        arad, 'arad')
        self.dlh = Util3d(model, (nlay, nrow, ncol), np.float32, dlh, 'dlh')
        self.dlv = Util3d(model, (nlay, nrow, ncol), np.float32, dlv, 'dlv')
        self.dth = Util3d(model, (nlay, nrow, ncol), np.float32, dth, 'dth')
        self.dtv = Util3d(model, (nlay, nrow, ncol), np.float32, dth, 'dtv')
        self.sconc = Util3d(model, (nlay, nrow, ncol), np.float32, sconc,
                            'sconc', )
        self.parent.add_package(self)
        return

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=extension, name=name,
                         unit_number=units, extra=extra, filenames=fname)

        nrow = model.nrow
        ncol = model.ncol
        nlay = model.nlay
        ncomp = model.ncomp
        mcomp = model.mcomp
        self.multiDiff = multiDiff
        self.al = Util3d(model, (nlay, nrow, ncol), np.float32, al, name='al',
                         locat=self.unit_number[0],
                         array_free_format=False)
        self.trpt = Util2d(model, (nlay,), np.float32, trpt, name='trpt',
                           locat=self.unit_number[0],
                           array_free_format=False)
        self.trpv = Util2d(model, (nlay,), np.float32, trpv, name='trpv',
                           locat=self.unit_number[0],
                           array_free_format=False)

        # Multi-species and multi-diffusion, hence the complexity
        self.dmcoef = []
        shape = (nlay, 1)
        utype = Util2d
        nmcomp = ncomp
        if multiDiff:
            shape = (nlay, nrow, ncol)

botm = np.linspace(top, botm, nlay)

        # Set values of all parameters
        self.heading = '# {} package for '.format(self.name[0]) + \
                       ' {}, '.format(model.version_types[model.version]) + \
                       'generated by Flopy.'
        self.laycbd = Util2d(model, (self.nlay,), np.int32, laycbd,
                             name='laycbd')
        self.laycbd[-1] = 0  # bottom layer must be zero
        self.delr = Util2d(model, (self.ncol,), np.float32, delr, name='delr',
                           locat=self.unit_number[0])
        self.delc = Util2d(model, (self.nrow,), np.float32, delc, name='delc',
                           locat=self.unit_number[0])
        self.top = Util2d(model, (self.nrow, self.ncol), np.float32,
                          top, name='model_top', locat=self.unit_number[0])
        self.botm = Util3d(model, (self.nlay + sum(self.laycbd),
                                   self.nrow, self.ncol), np.float32, botm,
                           'botm', locat=self.unit_number[0])
        self.perlen = Util2d(model, (self.nper,), np.float32, perlen,
                             name='perlen')
        self.nstp = Util2d(model, (self.nper,), np.int32, nstp, name='nstp')
        self.tsmult = Util2d(model, (self.nper,), np.float32, tsmult,
                             name='tsmult')
        self.steady = Util2d(model, (self.nper,), np.bool,
                             steady, name='steady')

        try:
            self.itmuni = int(itmuni)
        except:
            self.itmuni = ITMUNI[itmuni.lower()[0]]
        try:
            self.lenuni = int(lenuni)

#
        self.toeslope, self.tipslope, self.alpha, self.beta = toeslope, \
                                                              tipslope, \
                                                              alpha, \
                                                              beta
        self.nadptmx, self.nadptmn, self.adptfct = nadptmx, nadptmn, adptfct
        # Create arrays so that they have the correct size
        if self.istrat == 1:
            self.nu = Util2d(model, (self.nsrf + 1,), np.float32, nu,
                             name='nu')
        else:
            self.nu = Util2d(model, (self.nsrf + 2,), np.float32, nu,
                             name='nu')
        self.zeta = []
        for i in range(self.nsrf):
            self.zeta.append(Util3d(model, (nlay, nrow, ncol), np.float32,
                                    zeta[i], name='zeta_' + str(i + 1)))
        self.ssz = Util3d(model, (nlay, nrow, ncol), np.float32, ssz,
                          name='ssz')
        self.isource = Util3d(model, (nlay, nrow, ncol), np.int32, isource,
                              name='isource')
        #
        self.obsnam = obsnam
        self.obslrc = obslrc
        if nobs != 0:
            self.nobs = self.obslrc.shape[0]
        #
        self.parent.add_package(self)

self.wetfct = wetfct
        # Iteration interval for attempting to wet cells
        self.iwetit = iwetit
        # Flag that determines which equation is used to define the initial
        # head at cells that become wet
        self.ihdwet = ihdwet
        self.options = ' '
        if storagecoefficient:
            self.options = self.options + 'STORAGECOEFFICIENT '
        if constantcv: self.options = self.options + 'CONSTANTCV '
        if thickstrt: self.options = self.options + 'THICKSTRT '
        if nocvcorrection: self.options = self.options + 'NOCVCORRECTION '
        if novfc: self.options = self.options + 'NOVFC '
        self.hk = Util3d(model, (nlay, nrow, ncol), np.float32, hk, name='hk',
                         locat=self.unit_number[0])
        self.hani = Util3d(model, (nlay, nrow, ncol), np.float32, hani,
                           name='hani', locat=self.unit_number[0])
        keys = []
        for k in range(nlay):
            key = 'vka'
            if self.layvka[k] != 0:
                key = 'vani'
            keys.append(key)
        self.vka = Util3d(model, (nlay, nrow, ncol), np.float32, vka,
                          name=keys, locat=self.unit_number[0])
        tag = 'ss'
        if storagecoefficient:
            tag = 'storage'
        self.ss = Util3d(model, (nlay, nrow, ncol), np.float32, ss, name=tag,
                         locat=self.unit_number[0])
        self.sy = Util3d(model, (nlay, nrow, ncol), np.float32, sy, name='sy',
                         locat=self.unit_number[0])

# Create arrays so that they have the correct size
        if self.istrat == 1:
            self.nu = Util2d(model, (self.npln + 1,), np.float32, nu,
                             name='nu')
        else:
            self.nu = Util2d(model, (self.npln + 2,), np.float32, nu,
                             name='nu')
        self.zeta = []
        for i in range(nlay):
            # self.zeta.append( empty((nrow, ncol, self.npln)) )
            self.zeta
        for i in range(nlay):
            self.zeta.append(
                Util2d(model, (self.npln, nrow, ncol), np.float32, zeta[i],
                       name='zeta_' + str(i + 1)))
        self.ssz = Util3d(model, (nlay, nrow, ncol), np.float32, ssz,
                          name='ssz')
        self.isource = Util3d(model, (nlay, nrow, ncol), np.int, isource,
                              name='isource')
        self.parent.add_package(self)

name='Anisotropy factor', locat=self.unit_number[0])

        # item 1
        self.ipakcb = ipakcb
        self.hdry = hdry
        self.iwdflg = iwdflg
        self.wetfct = wetfct
        self.iwetit = iwetit
        self.ihdwet = ihdwet
        self.tran = Util3d(model, (nlay, nrow, ncol), np.float32, tran,
                           'Transmissivity', locat=self.unit_number[0])
        self.hy = Util3d(model, (nlay, nrow, ncol), np.float32, hy,
                         'Horizontal Hydraulic Conductivity',
                         locat=self.unit_number[0])
        if model.nlay > 1:
            self.vcont = Util3d(model, (nlay - 1, nrow, ncol), np.float32,
                                vcont,
                                'Vertical Conductance',
                                locat=self.unit_number[0])
        else:
            self.vcont = None
        self.sf1 = Util3d(model, (nlay, nrow, ncol), np.float32, sf1,
                          'Primary Storage Coefficient',
                          locat=self.unit_number[0])
        self.sf2 = Util3d(model, (nlay, nrow, ncol), np.float32, sf2,
                          'Secondary Storage Coefficient',
                          locat=self.unit_number[0])
        self.wetdry = Util3d(model, (nlay, nrow, ncol), np.float32, wetdry,
                             'WETDRY', locat=self.unit_number[0])
        self.parent.add_package(self)
        return

nrow, ncol, nlay, nper = self.parent.nrow_ncol_nlay_nper
        self.itrnsp = itrnsp
        self.ibctcb = ibctcb
        self.mcomp = mcomp
        self.ic_ibound_flg = ic_ibound_flg
        self.itvd = itvd
        self.iadsorb = iadsorb
        self.ict = ict
        self.cinact = cinact
        self.ciclose = ciclose
        self.idisp = idisp
        self.ixdisp = ixdisp
        self.diffnc = diffnc
        self.izod = izod
        self.ifod = ifod
        self.icbund = Util3d(model, (nlay, nrow, ncol), np.float32, icbund,
                             'icbund', )
        self.porosity = Util3d(model, (nlay, nrow, ncol), np.float32,
                               porosity, 'porosity')
        # self.arad = Util2d(model, (1, nja), np.float32,
        #                        arad, 'arad')
        self.dlh = Util3d(model, (nlay, nrow, ncol), np.float32, dlh, 'dlh')
        self.dlv = Util3d(model, (nlay, nrow, ncol), np.float32, dlv, 'dlv')
        self.dth = Util3d(model, (nlay, nrow, ncol), np.float32, dth, 'dth')
        self.dtv = Util3d(model, (nlay, nrow, ncol), np.float32, dth, 'dtv')
        self.sconc = Util3d(model, (nlay, nrow, ncol), np.float32, sconc,
                            'sconc', )
        self.parent.add_package(self)
        return

def thickness(self):
        """
        Get a Util3d array of cell thicknesses.

        Returns
        -------
        thickness : util3d array of floats (nlay, nrow, ncol)

        """
        # return self.__thickness
        thk = []
        thk.append(self.top - self.botm[0])
        for k in range(1, self.nlay + sum(self.laycbd)):
            thk.append(self.botm[k - 1] - self.botm[k])
        return Util3d(self.parent, (self.nlay + sum(self.laycbd),
                                    self.nrow, self.ncol), np.float32,
                      thk, name='thickness')