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

sp1_wells.loc[:, "flux"] = 2200.0
    sp1_wells = sp1_wells.loc[:, ["k", "i", "j", "flux"]].values.tolist()

    sp2_wells = sp1_wells.copy()
    sp2_wells.append([1, 8, 9, -72000.0])
    sp2_wells.append([3, 11, 6, -72000.0])

    flopy.modflow.ModflowWel(ml, stress_period_data=
                            {0: sp1_wells, 1: sp2_wells, 2: sp1_wells})

    flopy.modflow.ModflowSwt(ml,iswtoc=1,nsystm=4,sgs=2.0,sgm=1.7,
                               lnwt=[0,1,2,3],thick=[45,70,50,90],icrcc=0,
                               cr=0.01,cc=0.25,istpcs=1,pcsoff=15.0,
                               void=0.82,ithk=1,ivoid=1)
    flopy.modflow.ModflowOc(ml,stress_period_data={(0,0):["save head","save budget"]})
    flopy.modflow.ModflowPcg(ml,hclose=0.01,rclose=1.0)

    ml.write_input()

    ml1 = flopy.modflow.Modflow.load(ml.namefile,model_ws=model_ws,
                                     verbose=True,exe_name="mf2005",load_only=["SWT"],forgive=False)

    return ml

# creating MODFLOW model
    
    model_ws = os.path.join('.', 'temp', 't068a')
    modelname = 'model_mf'
    
    mf = flopy.modflow.Modflow(modelname, model_ws=model_ws,
                               exe_name=mf_exe_name)
    dis = flopy.modflow.ModflowDis(mf, nlay=nlay, nrow=nrow, ncol=ncol, 
                                   perlen=perlen, nper=nper, botm=botm, top=top, 
                                   steady=False)
    
    bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=top)
    lpf = flopy.modflow.ModflowLpf(mf, hk=100, vka=100, ss=0.00001, sy=0.1)
    oc = flopy.modflow.ModflowOc(mf)
    pcg = flopy.modflow.ModflowPcg(mf)
    lmt = flopy.modflow.ModflowLmt(mf)
    
    # recharge
    rchrate = {}
    rchrate[0] = 0.0
    rchrate[5] = 0.001
    rchrate[6] = 0.0
    
    rch = flopy.modflow.ModflowRch(mf, rech=rchrate, nrchop=3)
    
    # define itype
    itype = flopy.mt3d.Mt3dSsm.itype_dict()
    ssm_data = {}
    
    # Now for the point of this test: Enter SSM data sometime
    # after the first stress period (this was crashing flopy

flopy.modflow.ModflowBas(m, ibound=1, strt=top)
    # recharge
    flopy.modflow.ModflowRch(m, ipakcb=iu_cbc, rech=rch, nrchop=1)
    # wel
    wd = [i for i in wel_loc] + [wel_q]
    flopy.modflow.ModflowWel(m, ipakcb=iu_cbc, stress_period_data={0: wd})
    # river
    rd = []
    for i in range(nrow):
        rd.append([0, i, ncol - 1, riv_h, riv_c, riv_z])
    flopy.modflow.ModflowRiv(m, ipakcb=iu_cbc, stress_period_data={0: rd})
    # output control
    flopy.modflow.ModflowOc(m, stress_period_data={(0, 0): ['save head',
                                                            'save budget',
                                                            'print head']})
    flopy.modflow.ModflowPcg(m, hclose=1e-6, rclose=1e-3, iter1=100, mxiter=50)

    m.write_input()

    if run:
        success, buff = m.run_model()
        assert success, 'mf2005 model did not run'

    # create modpath files
    exe_name = exe_names['mp7']
    mp = flopy.modpath.Modpath7(modelname=nm + '_mp', flowmodel=m,
                                exe_name=exe_name, model_ws=ws)
    mpbas = flopy.modpath.Modpath7Bas(mp, porosity=0.1,
                                      defaultiface=defaultiface)
    mpsim = flopy.modpath.Modpath7Sim(mp, simulationtype='combined',
                                      trackingdirection='forward',
                                      weaksinkoption='pass_through',

model_ws=workspace)
    # create flopy modflow package objects
    discret = flopy.modflow.ModflowDis(ml, nlay=nlay, nrow=nrow, ncol=ncol,
                                       delr=delr, delc=delc,
                                       top=0, botm=[-40.0],
                                       perlen=400, nstp=200)
    bas = flopy.modflow.ModflowBas(ml, ibound=ibound, strt=0.0)
    lpf = flopy.modflow.ModflowLpf(ml, hk=2., vka=2.0, vkcb=0, laytyp=0,
                                   layavg=0)
    wel = flopy.modflow.ModflowWel(ml, stress_period_data={0: [(0, 0, 0, 1)]})
    swi = flopy.modflow.ModflowSwi2(ml, iswizt=55, npln=1, istrat=1,
                                    toeslope=0.2, tipslope=0.2, nu=[0, 0.025],
                                    zeta=z, ssz=0.2, isource=isource,
                                    nsolver=1)
    oc = flopy.modflow.ModflowOc(ml, stress_period_data=ocdict)
    pcg = flopy.modflow.ModflowPcg(ml)
    # create model files
    ml.write_input()
    # run the model
    m = ml.run_model(silent=False)
    # read model heads
    headfile = os.path.join(workspace, '{}.hds'.format(modelname))
    hdobj = flopy.utils.HeadFile(headfile)
    head = hdobj.get_alldata()
    head = np.array(head)
    # read model zeta
    zetafile = os.path.join(workspace, '{}.zta'.format(modelname))
    zobj = flopy.utils.CellBudgetFile(zetafile)
    zkstpkper = zobj.get_kstpkper()
    zeta = []
    for kk in zkstpkper:
        zeta.append(zobj.get_data(kstpkper=kk, text='      ZETASRF  1')[0])

# Add DIS package to the MODFLOW model
    dis = flopy.modflow.ModflowDis(m, nlay, nrow, ncol, nper=1, delr=delr,
                                   delc=delc, laycbd=0, top=henry_top,
                                   botm=henry_botm, perlen=1.5, nstp=15)

    # Variables for the BAS package
    ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)
    ibound[:, :, -1] = -1
    bas = flopy.modflow.ModflowBas(m, ibound, 0)

    # Add LPF package to the MODFLOW model
    lpf = flopy.modflow.ModflowLpf(m, hk=hk, vka=hk, ipakcb=53)

    # Add PCG Package to the MODFLOW model
    pcg = flopy.modflow.ModflowPcg(m, hclose=1.e-8)

    # Add OC package to the MODFLOW model
    oc = flopy.modflow.ModflowOc(m,
                                 stress_period_data={
                                     (0, 0): ['save head', 'save budget']},
                                 compact=True)

    # Create WEL and SSM data
    itype = flopy.mt3d.Mt3dSsm.itype_dict()
    wel_data = {}
    ssm_data = {}
    wel_sp1 = []
    ssm_sp1 = []
    for k in range(nlay):
        wel_sp1.append([k, 0, 0, qinflow / nlay])
        ssm_sp1.append([k, 0, 0, 0., itype['WEL']])

nper = 3
perlen = [1, 100, 100]
nstp = [1, 100, 100]
steady = [True, False, False]

# Flopy objects
modelname = 'tutorial2'
mf = flopy.modflow.Modflow(modelname, exe_name='mf2005')
dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc,
                               top=ztop, botm=botm[1:],
                               nper=nper, perlen=perlen, nstp=nstp,
                               steady=steady)
bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)
lpf = flopy.modflow.ModflowLpf(mf, hk=hk, vka=vka, sy=sy, ss=ss, laytyp=laytyp,
                               ipakcb=53)
pcg = flopy.modflow.ModflowPcg(mf)

# Make list for stress period 1
stageleft = 10.
stageright = 10.
bound_sp1 = []
for il in range(nlay):
    condleft = hk * (stageleft - zbot) * delc
    condright = hk * (stageright - zbot) * delc
    for ir in range(nrow):
        bound_sp1.append([il, ir, 0, stageleft, condleft])
        bound_sp1.append([il, ir, ncol - 1, stageright, condright])
print('Adding ', len(bound_sp1), 'GHBs for stress period 1.')

# Make list for stress period 2
stageleft = 10.
stageright = 0.

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')

"wel": flopy.modflow.ModflowWel,
            "mnw1": flopy.modflow.ModflowMnw1,
            "mnw2": flopy.modflow.ModflowMnw2,
            "mnwi": flopy.modflow.ModflowMnwi,
            "drn": flopy.modflow.ModflowDrn,
            "drt": flopy.modflow.ModflowDrt,
            "rch": flopy.modflow.ModflowRch,
            "evt": flopy.modflow.ModflowEvt,
            "ghb": flopy.modflow.ModflowGhb,
            "gmg": flopy.modflow.ModflowGmg,
            "lmt6": flopy.modflow.ModflowLmt,
            "lmt7": flopy.modflow.ModflowLmt,
            "riv": flopy.modflow.ModflowRiv,
            "str": flopy.modflow.ModflowStr,
            "swi2": flopy.modflow.ModflowSwi2,
            "pcg": flopy.modflow.ModflowPcg,
            "pcgn": flopy.modflow.ModflowPcgn,
            "nwt": flopy.modflow.ModflowNwt,
            "pks": flopy.modflow.ModflowPks,
            "sms": flopy.modflow.ModflowSms,
            "sfr": flopy.modflow.ModflowSfr2,
            "lak": flopy.modflow.ModflowLak,
            "gage": flopy.modflow.ModflowGage,
            "sip": flopy.modflow.ModflowSip,
            "sor": flopy.modflow.ModflowSor,
            "de4": flopy.modflow.ModflowDe4,
            "oc": flopy.modflow.ModflowOc,
            "uzf": flopy.modflow.ModflowUzf1,
            "upw": flopy.modflow.ModflowUpw,
            "sub": flopy.modflow.ModflowSub,
            "swt": flopy.modflow.ModflowSwt,
            "hyd": flopy.modflow.ModflowHyd,