How to use the flopy.modflow.ModflowOc function in flopy

# control file
    # Add OC package to the MODFLOW model
    options = ['PRINT HEAD', 'PRINT DRAWDOWN', 'PRINT BUDGET',
               'SAVE HEAD', 'SAVE DRAWDOWN', 'SAVE BUDGET',
               'SAVE IBOUND', 'DDREFERENCE']
    idx = 0
    spd = dict()
    for sp in mf.dis.nstp.array:
        stress_period = idx
        step = sp - 1
        ke = (stress_period, step)
        idx = idx + 1
        spd[ke] = [options[3], options[2], options[5]]
    mf.remove_package('OC')
    oc = flopy.modflow.ModflowOc(mf, stress_period_data=spd, cboufm='(20i5)')
    oc.write_file()

    # write hob files
    hob_df = pd.read_csv(r"misc_hob.csv", header=None)
    start_time = hob_df.values[0][0]
    end_time = hob_df.values[1][0]

    #hobs

    mf_name2 = r".\mf_EnS\pp_model.nam"
    mf2 = flopy.modflow.Modflow.load(mf_name2, load_only=['DIS', 'BAS6','HOB'])

    obs_data = mf2.hob.obs_data
    if len(steady) > 1:
        time_index = np.arange(start_time+1,end_time+1)
    else:

nam_file = "freyberg.nam"
    mo = flopy.modflow.Modflow.load(nam_file, model_ws=org_model_ws, check=False,forgive=False)

    perlen = np.ones((365))
    m = flopy.modflow.Modflow("freyberg_transient",model_ws=os.path.join("da","freyberg","truth"),version="mfnwt",
                              external_path=".")

    flopy.modflow.ModflowDis(m,nrow=mo.nrow,ncol=mo.ncol,nlay=1,delr=mo.dis.delr,
                             delc=mo.dis.delc,top=mo.dis.top,botm=mo.dis.botm[-1],nper=len(perlen),perlen=perlen)
    flopy.modflow.ModflowBas(m,ibound=mo.bas6.ibound[0],strt=mo.bas6.strt[0])
    flopy.modflow.ModflowUpw(m,laytyp=mo.upw.laytyp,hk=mo.upw.hk[0],vka=mo.upw.vka[0],ss=0.00001,sy=0.01)
    flopy.modflow.ModflowNwt(m)
    oc_data = {}
    for iper in range(m.nper):
        oc_data[iper,0] = ["save head","save budget"]
    flopy.modflow.ModflowOc(m,stress_period_data=oc_data)
    flopy.modflow.ModflowRch(m,rech=mo.rch.rech.array[0])
    wel_data = mo.wel.stress_period_data[0]
    wel_data["k"][:] = 0
    flopy.modflow.ModflowWel(m,stress_period_data={0:wel_data})
    flopy.modflow.ModflowSfr2(m,nstrm=mo.sfr.nstrm,nss=mo.sfr.nss,istcb2=90,segment_data=mo.sfr.segment_data,reach_data=mo.sfr.reach_data)
    m.write_input()
    pyemu.os_utils.run("mfnwt {0}.nam".format(m.name),cwd=m.model_ws)
    return m

file = '{}_ib{:02d}.ref'.format(name, kdx + 1)
    files.append(file)
    hfiles.append(hfile)
    np.savetxt(file, ibound[kdx, :, :], fmt='%5d')

bas = flopy.modflow.ModflowBas(ml, ibound=files, strt=hfiles)

# The aquifer properties (really only the hydraulic conductivity) are defined with the
# LPF package.
lpf = flopy.modflow.ModflowLpf(ml, hk=Kh)

# Finally, we need to specify the solver we want to use (PCG with default values), and the
# output control (using the default values). Then we are ready to write all MODFLOW input
# files and run MODFLOW.
pcg = flopy.modflow.ModflowPcg(ml)
oc = flopy.modflow.ModflowOc(ml)
ml.write_input()
ml.run_model()

# change back to the starting directory
os.chdir(cwdpth)

# Once the model has terminated normally, we can read the heads file. First, a link to the heads
# file is created with `HeadFile`. The link can then be accessed with the `get_data` function, by
# specifying, in this case, the step number and period number for which we want to retrieve data.
# A three-dimensional array is returned of size `nlay, nrow, ncol`. Matplotlib contouring functions
# are used to make contours of the layers or a cross-section.
hds = flopy.utils.HeadFile(os.path.join(workspace, name + '.hds'))
h = hds.get_data(kstpkper=(0, 0))
x = y = np.linspace(0, L, N)
c = plt.contour(x, y, h[0], np.arange(90, 100.1, 0.2))
plt.clabel(c, fmt='%2.1f')

wel_sp1 = [[0, nrow/2 - 1, ncol/2 - 1, 0.]]
wel_sp2 = [[0, nrow/2 - 1, ncol/2 - 1, 0.]]
wel_sp3 = [[0, nrow/2 - 1, ncol/2 - 1, pumping_rate]]
stress_period_data = {0: wel_sp1, 1: wel_sp2, 2: wel_sp3}
wel = flopy.modflow.ModflowWel(mf, stress_period_data=stress_period_data)

# Output control
stress_period_data = {}
for kper in range(nper):
    for kstp in range(nstp[kper]):
        stress_period_data[(kper, kstp)] = ['save head',
                                            'save drawdown',
                                            'save budget',
                                            'print head',
                                            'print budget']
oc = flopy.modflow.ModflowOc(mf, stress_period_data=stress_period_data,
                             compact=True)

# Write the model input files
mf.write_input()

# Run the model
success, mfoutput = mf.run_model(silent=False, pause=False)
if not success:
    raise Exception('MODFLOW did not terminate normally.')


# Imports
import matplotlib.pyplot as plt
import flopy.utils.binaryfile as bf

# Create the headfile and budget file objects

print('creating...', modelname)
    ml = flopy.modflow.Modflow(modelname, version='mf2005', exe_name=mf_name,
                               model_ws=dirs[0])
    discret = flopy.modflow.ModflowDis(ml, nlay=1, ncol=ncol, nrow=nrow,
                                       delr=delr,
                                       delc=1,
                                       top=0, botm=[-40.0],
                                       nper=nper, perlen=perlen, nstp=nstp)
    bas = flopy.modflow.ModflowBas(ml, ibound=ibound, strt=0.05)
    bcf = flopy.modflow.ModflowBcf(ml, laycon=0, tran=2 * 40)
    swi = flopy.modflow.ModflowSwi2(ml, iswizt=55, nsrf=nsurf, istrat=1,
                                    toeslope=0.2, tipslope=0.2,
                                    nu=[0, 0.0125, 0.025],
                                    zeta=z, ssz=ssz, isource=isource,
                                    nsolver=1)
    oc = flopy.modflow.ModflowOc(ml,
                                 stress_period_data={(0, 999): ['save head']})
    pcg = flopy.modflow.ModflowPcg(ml)
    ml.write_input()
    # run stratified model
    if not skipRuns:
        m = ml.run_model(silent=False)
    # read stratified results
    zetafile = os.path.join(dirs[0], '{}.zta'.format(modelname))
    zobj = flopy.utils.CellBudgetFile(zetafile)
    zkstpkper = zobj.get_kstpkper()
    zeta = zobj.get_data(kstpkper=zkstpkper[-1], text='ZETASRF  1')[0]
    zeta2 = zobj.get_data(kstpkper=zkstpkper[-1], text='ZETASRF  2')[0]
    #
    # vd model
    modelname = 'swiex2_vd'
    print('creating...', modelname)

bot = np.linspace(-H/Nlay,-H,Nlay)
delrow = delcol = L/(N-1)
dis = fmf.ModflowDis(ml,nlay=Nlay,nrow=N,ncol=N,delr=delrow,delc=delcol,top=0.0,botm=bot,laycbd=0)

Nhalf = (N-1)/2
ibound = np.ones((Nlay,N,N))
ibound[:,0,:] = -1; ibound[:,-1,:] = -1; ibound[:,:,0] = -1; ibound[:,:,-1] = -1
ibound[0,Nhalf,Nhalf] = -1
start = h1 * np.ones((N,N))
start[Nhalf,Nhalf] = h2
bas = fmf.ModflowBas(ml,ibound=ibound,strt=start)
print "o"
lpf = fmf.ModflowLpf(ml, hk=k)
pcg = fmf.ModflowPcg(ml)
oc = fmf.ModflowOc(ml)
ml.write_input()
ml.run_model()

head_file = '/home/kiruba/PycharmProjects/area_of_curve/hydrology/hydrology/mf_files/lake_example.hds'
hds = fut.HeadFile(head_file)
h = hds.get_data(kstpkper=(1,1))
print len(h)
x = y = np.linspace(0, L, N)
c = plt.contour(x, y, h[0], np.arange(90,100.1,0.2))
plt.clabel(c, fmt='%2.1f')
plt.axis('scaled')
plt.show()
c = plt.contour(x,y,h[-1],np.arange(90,100.1,0.2))
plt.clabel(c,fmt='%1.1f')
plt.axis('scaled')
plt.show()

(1,30): [],
                                  (1,39): ['save head', 'save budget', 'print budget'], 
                                  (1,40): [],
                                  (1,49): ['save head', 'save budget', 'print budget'], 
                                  (1,50): [],
                                  (1,59): ['save head', 'save budget', 'print budget'], 
                                  (1,60): [],
                                  (1,69): ['save head', 'save budget', 'print budget'], 
                                  (1,70): [],
                                  (1,79): ['save head', 'save budget', 'print budget'], 
                                  (1,80): [],
                                  (1,89): ['save head', 'save budget', 'print budget'], 
                                  (1,90): [],
                                  (1,99): ['save head', 'save budget', 'print budget'], 
                                  (1,100): []}
            oc = flopy.modflow.ModflowOc(ml, stress_period_data=stress_period_data)
            # write the input files
            ml.write_input()
        else:
            if not os.path.exists(cf_pths[idx]):
                os.makedirs(cf_pths[idx])
            filelist = [f for f in os.listdir(cf_pths[0])]
            sys.stdout.write('copying files from {} to {}\n'.format(cf_pths[0], cf_pths[idx]))
            for f in filelist:
                if os.path.splitext(f)[1].lower() in exclude:
                    continue
                src = os.path.join(cf_pths[0], f)
                dst = os.path.join(cf_pths[idx], f)
                shutil.copyfile(src, dst)

    return ml, cf_pths