In [13]:
import os
import numpy as np
import flopy
Setup the basic dir structure. Since we are using external arrays, we have to change into the model_ws dir.
In [14]:
starting_ws = os.getcwd()
model_ws = "henry"
if not os.path.exists(model_ws):
os.mkdir(model_ws)
# change to base_dir to write the model files
os.chdir(model_ws)
In [15]:
ml = flopy.modflow.Modflow(modelname="henry", external_path="ref")
ml.array_free_format = False
DIS package
In [16]:
nlay, nrow, ncol = 20, 1, 120
delr, delc = 0.1, 1.0
delz = 0.1
botm = np.linspace(-delz, -float(nlay)*delz, nlay)
dis = flopy.modflow.ModflowDis(ml, nrow=nrow, ncol=ncol,
nlay=nlay, nper=2, perlen=[1,0.33],
steady=[True, False], delr=delr,
delc=delc, botm=botm, top=0.0)
BAS package
In [17]:
bas = flopy.modflow.ModflowBas(ml, ibound=1, strt=0.0)
LPF package using external arrays for hk and vka
In [18]:
hk = np.zeros((nlay, nrow, ncol))+200.0
lpf = flopy.modflow.ModflowLpf(ml, hk=hk, vka=hk, ss=1.0E-5, sy=0.15, laytyp=1)
WEL, GHB and PCG packages. WEL's are set along the left side of the domain to represent freshwater inflow and GHB's are set along the right side of the domain to represent a saltwater boundary. Also track SSM entries for later
In [19]:
well_data, ghb_data, ssm_data = {}, {}, {}
wel_1, ghb_1, ssm_1 = [], [], []
itype = flopy.mt3d.Mt3dSsm.itype_dict()
for k in xrange(nlay):
wel_1.append([k, 0, 0, 0.2])
ssm_1.append([k, 0, 0, 0.0, itype["WEL"]])
ghb_1.append([k, 0, ncol-1, 0.0, 1.0E+3])
ssm_1.append([k, 0, ncol-1, 1.0, itype["GHB"]])
ghb_data[0] = ghb_1
well_data[0] = wel_1
# turn off the freshwater flux in sp 2
well_data[1] = 0
ssm_data[0] = ssm_1
ghb = flopy.modflow.ModflowGhb(ml, stress_period_data=ghb_data)
wel = flopy.modflow.ModflowWel(ml, stress_period_data=well_data)
pcg = flopy.modflow.ModflowPcg(ml,)
write the modflow input
In [20]:
ml.write_input()
In [21]:
mt = flopy.mt3d.Mt3dms("henry", namefile_ext="nam_mt", modflowmodel=ml,
external_path="ref")#,model_ws=model_ws)
btn = flopy.mt3d.Mt3dBtn(mt)
adv = flopy.mt3d.Mt3dAdv(mt, mixelm=0)
ssm = flopy.mt3d.Mt3dSsm(mt, stress_period_data=ssm_data)
gcg = flopy.mt3d.Mt3dGcg(mt)
mt.write_input()
swt = flopy.seawat.Seawat(modelname="henry", namefile_ext="nam_swt",
modflowmodel=ml, mt3dmsmodel=mt)#,model_ws=model_ws)
vdf = flopy.seawat.SeawatVdf(swt, densemin=0.0, densemax=0.0, denseslp=24.5)
swt.write_input()
change back to the original dir
In [22]:
os.chdir(starting_ws)
verify that we can load the modflow
In [23]:
ml_loaded = flopy.modflow.Modflow.load("henry.nam", model_ws=model_ws)
ml_loaded.change_model_ws("new_henry")
ml_loaded.write_input()
In [24]:
print ml_loaded.get_package("lpf").hk.array