<s> from __future__ import print_function <EOL> import os <EOL> import sys <EOL> import math <EOL> import numpy as np <EOL> import flopy . modflow as mf <EOL> import flopy . mt3d as mt3 <EOL> import flopy . seawat as swt <EOL> import flopy . utils as fu <EOL> import matplotlib . pyplot as plt <EOL> updates = { '<STR_LIT>' : [ '<STR_LIT>' , '<STR_LIT>' ] , <EOL> '<STR_LIT>' : '<STR_LIT>' , <EOL> '<STR_LIT>' : <NUM_LIT:0> , <EOL> '<STR_LIT>' : <NUM_LIT> , <EOL> '<STR_LIT>' : <NUM_LIT:7> , <EOL> '<STR_LIT>' : <NUM_LIT:8> , <EOL> '<STR_LIT>' : <NUM_LIT:7> , <EOL> '<STR_LIT>' : <NUM_LIT:7> } <EOL> plt . rcParams . update ( updates ) <EOL> cleanFiles = False <EOL> skipRuns = False <EOL> fext = '<STR_LIT>' <EOL> narg = len ( sys . argv ) <EOL> iarg = <NUM_LIT:0> <EOL> if narg > <NUM_LIT:1> : <EOL> while iarg < narg - <NUM_LIT:1> : <EOL> iarg += <NUM_LIT:1> <EOL> basearg = sys . argv [ iarg ] . lower ( ) <EOL> if basearg == '<STR_LIT>' : <EOL> cleanFiles = True <EOL> elif basearg == '<STR_LIT>' : <EOL> skipRuns = True <EOL> elif basearg == '<STR_LIT>' : <EOL> fext = '<STR_LIT>' <EOL> dirs = [ os . path . join ( '<STR_LIT>' ) , os . path . join ( '<STR_LIT>' ) ] <EOL> if cleanFiles : <EOL> print ( '<STR_LIT>' ) <EOL> print ( '<STR_LIT>' ) <EOL> file_dict = { } <EOL> file_dict [ '<STR_LIT:.>' ] = os . listdir ( '<STR_LIT:.>' ) <EOL> file_dict [ dirs [ <NUM_LIT:0> ] ] = os . listdir ( dirs [ <NUM_LIT:0> ] ) <EOL> file_dict [ dirs [ <NUM_LIT:1> ] ] = os . listdir ( dirs [ <NUM_LIT:1> ] ) <EOL> for key , files in file_dict . items ( ) : <EOL> for f in files : <EOL> if os . path . isdir ( f ) : <EOL> continue <EOL> if '<STR_LIT>' != os . path . splitext ( f ) [ <NUM_LIT:1> ] . lower ( ) : <EOL> print ( '<STR_LIT>' . format ( os . path . basename ( f ) ) ) <EOL> os . remove ( os . path . join ( key , f ) ) <EOL> for d in dirs : <EOL> if os . path . exists ( d ) : <EOL> os . rmdir ( d ) <EOL> sys . exit ( <NUM_LIT:1> ) <EOL> for d in dirs : <EOL> if not os . path . exists ( d ) : <EOL> os . mkdir ( d ) <EOL> nlay = <NUM_LIT:6> <EOL> nrow = <NUM_LIT:1> <EOL> ncol = <NUM_LIT> <EOL> delr = np . zeros ( ( ncol ) , np . float ) <EOL> delc = <NUM_LIT:1.> <EOL> r = np . zeros ( ( ncol ) , np . float ) <EOL> x = np . zeros ( ( ncol ) , np . float ) <EOL> edge = np . zeros ( ( ncol ) , np . float ) <EOL> dx = <NUM_LIT> <EOL> for i in range ( <NUM_LIT:0> , ncol ) : <EOL> delr [ i ] = dx <EOL> r [ <NUM_LIT:0> ] = delr [ <NUM_LIT:0> ] / <NUM_LIT> <EOL> for i in range ( <NUM_LIT:1> , ncol ) : <EOL> r [ i ] = r [ i - <NUM_LIT:1> ] + ( delr [ i - <NUM_LIT:1> ] + delr [ i ] ) / <NUM_LIT> <EOL> x [ <NUM_LIT:0> ] = delr [ <NUM_LIT:0> ] / <NUM_LIT> <EOL> for i in range ( <NUM_LIT:1> , ncol ) : <EOL> x [ i ] = x [ i - <NUM_LIT:1> ] + ( delr [ i - <NUM_LIT:1> ] + delr [ i ] ) / <NUM_LIT> <EOL> edge [ <NUM_LIT:0> ] = delr [ <NUM_LIT:0> ] <EOL> for i in range ( <NUM_LIT:1> , ncol ) : <EOL> edge [ i ] = edge [ i - <NUM_LIT:1> ] + delr [ i ] <EOL> nper = <NUM_LIT:2> <EOL> perlen = [ <NUM_LIT> , <NUM_LIT> ] <EOL> nstp = [ <NUM_LIT> , <NUM_LIT> ] <EOL> steady = True <EOL> nsave_zeta = <NUM_LIT:8> <EOL> ndecay = <NUM_LIT:4> <EOL> ibound = np . ones ( ( nlay , nrow , ncol ) , np . int ) <EOL> for k in range ( <NUM_LIT:0> , nlay ) : <EOL> ibound [ k , <NUM_LIT:0> , ncol - <NUM_LIT:1> ] = - <NUM_LIT:1> <EOL> bot = np . zeros ( ( nlay , nrow , ncol ) , np . float ) <EOL> dz = <NUM_LIT> / float ( nlay - <NUM_LIT:1> ) <EOL> zall = - np . arange ( <NUM_LIT:0> , <NUM_LIT:100> + dz , dz ) <EOL> zall = np . append ( zall , - <NUM_LIT> ) <EOL> tb = - np . arange ( dz , <NUM_LIT:100> + dz , dz ) <EOL> tb = np . append ( tb , - <NUM_LIT> ) <EOL> for k in range ( <NUM_LIT:0> , nlay ) : <EOL> for i in range ( <NUM_LIT:0> , ncol ) : <EOL> bot [ k , <NUM_LIT:0> , i ] = tb [ k ] <EOL> isource = np . zeros ( ( nlay , nrow , ncol ) , np . int ) <EOL> isource [ : , <NUM_LIT:0> , ncol - <NUM_LIT:1> ] = <NUM_LIT:1> <EOL> isource [ nlay - <NUM_LIT:1> , <NUM_LIT:0> , ncol - <NUM_LIT:1> ] = <NUM_LIT:2> <EOL> khb = ( <NUM_LIT> * <NUM_LIT> * <NUM_LIT> / <NUM_LIT> ) * <NUM_LIT> * <NUM_LIT> * <NUM_LIT> <EOL> kvb = ( <NUM_LIT> * <NUM_LIT> * <NUM_LIT> / <NUM_LIT> ) * <NUM_LIT> * <NUM_LIT> * <NUM_LIT> <EOL> ssb = <NUM_LIT> <EOL> sszb = <NUM_LIT> <EOL> kh = np . zeros ( ( nlay , nrow , ncol ) , np . float ) <EOL> kv = np . zeros ( ( nlay , nrow , ncol ) , np . float ) <EOL> ss = np . zeros ( ( nlay , nrow , ncol ) , np . float ) <EOL> ssz = np . zeros ( ( nlay , nrow , ncol ) , np . float ) <EOL> for k in range ( <NUM_LIT:0> , nlay ) : <EOL> for i in range ( <NUM_LIT:0> , ncol ) : <EOL> f = r [ i ] * <NUM_LIT> * math . pi <EOL> kh [ k , <NUM_LIT:0> , i ] = khb * f <EOL> kv [ k , <NUM_LIT:0> , i ] = kvb * f <EOL> ss [ k , <NUM_LIT:0> , i ] = ssb * f <EOL> ssz [ k , <NUM_LIT:0> , i ] = sszb * f <EOL> z = np . ones ( ( nlay ) , np . float ) <EOL> z = - <NUM_LIT> * z <EOL> nwell = <NUM_LIT:1> <EOL> for k in range ( <NUM_LIT:0> , nlay ) : <EOL> if zall [ k ] > - <NUM_LIT> and zall [ k + <NUM_LIT:1> ] <= - <NUM_LIT:20> : <EOL> nwell = k + <NUM_LIT:1> <EOL> print ( '<STR_LIT>' . format ( nlay , dz , nwell ) ) <EOL> wellQ = - <NUM_LIT> <EOL> wellbtm = - <NUM_LIT> <EOL> wellQpm = wellQ / abs ( wellbtm ) <EOL> well_data = { } <EOL> for ip in range ( <NUM_LIT:0> , nper ) : <EOL> welllist = np . zeros ( ( nwell , <NUM_LIT:4> ) , np . float ) <EOL> for iw in range ( <NUM_LIT:0> , nwell ) : <EOL> if ip == <NUM_LIT:0> : <EOL> b = zall [ iw ] - zall [ iw + <NUM_LIT:1> ] <EOL> if zall [ iw + <NUM_LIT:1> ] < wellbtm : <EOL> b = zall [ iw ] - wellbtm <EOL> q = wellQpm * b <EOL> else : <EOL> q = <NUM_LIT:0.0> <EOL> welllist [ iw , <NUM_LIT:0> ] = iw <EOL> welllist [ iw , <NUM_LIT:1> ] = <NUM_LIT:0> <EOL> welllist [ iw , <NUM_LIT:2> ] = <NUM_LIT:0> <EOL> welllist [ iw , <NUM_LIT:3> ] = q <EOL> well_data [ ip ] = welllist . copy ( ) <EOL> ihead = np . zeros ( ( nlay ) , np . float ) <EOL> savewords = [ ] <EOL> for i in range ( <NUM_LIT:0> , nper ) : <EOL> icnt = <NUM_LIT:0> <EOL> for j in range ( <NUM_LIT:0> , nstp [ i ] ) : <EOL> icnt += <NUM_LIT:1> <EOL> savebudget = False <EOL> savehead = False <EOL> if icnt == <NUM_LIT> : <EOL> savebudget = True <EOL> savehead = True <EOL> icnt = <NUM_LIT:0> <EOL> if j < <NUM_LIT:10> : <EOL> savebudget = True <EOL> if savebudget == True or savehead == True : <EOL> twords = [ i + <NUM_LIT:1> , j + <NUM_LIT:1> ] <EOL> if savebudget == True : <EOL> twords . append ( '<STR_LIT>' ) <EOL> if savehead == True : <EOL> twords . append ( '<STR_LIT>' ) <EOL> savewords . append ( twords ) <EOL> solver2params = { '<STR_LIT>' : <NUM_LIT:100> , '<STR_LIT>' : <NUM_LIT:20> , '<STR_LIT>' : <NUM_LIT:1> , '<STR_LIT>' : <NUM_LIT> , '<STR_LIT>' : <NUM_LIT> , '<STR_LIT>' : <NUM_LIT:1.0> , <EOL> '<STR_LIT>' : <NUM_LIT:2> , '<STR_LIT>' : <NUM_LIT:1.0> , '<STR_LIT>' : <NUM_LIT:1.0> } <EOL> modelname = '<STR_LIT>' <EOL> mf_name = '<STR_LIT>' <EOL> if not skipRuns : <EOL> ml = mf . Modflow ( modelname , version = '<STR_LIT>' , exe_name = mf_name , model_ws = dirs [ <NUM_LIT:0> ] ) <EOL> discret = mf . ModflowDis ( ml , nrow = nrow , ncol = ncol , nlay = nlay , delr = delr , delc = delc , <EOL> top = <NUM_LIT:0> , botm = bot , laycbd = <NUM_LIT:0> , nper = nper , <EOL> perlen = perlen , nstp = nstp , steady = steady ) <EOL> bas = mf . ModflowBas ( ml , ibound = ibound , strt = ihead ) <EOL> lpf = mf . ModflowLpf ( ml , hk = kh , vka = kv , ss = ss , sy = ssz , vkcb = <NUM_LIT:0> , laytyp = <NUM_LIT:0> , layavg = <NUM_LIT:1> ) <EOL> wel = mf . ModflowWel ( ml , stress_period_data = well_data ) <EOL> swi = mf . ModflowSwi2 ( ml , npln = <NUM_LIT:1> , istrat = <NUM_LIT:1> , toeslope = <NUM_LIT> , tipslope = <NUM_LIT> , nu = [ <NUM_LIT:0> , <NUM_LIT> ] , zeta = z , ssz = ssz , isource = isource , nsolver = <NUM_LIT:2> , solver2params = solver2params ) <EOL> oc = mf . ModflowOc88 ( ml , words = savewords ) <EOL> pcg = mf . ModflowPcg ( ml , hclose = <NUM_LIT> , rclose = <NUM_LIT> , mxiter = <NUM_LIT:100> , iter1 = <NUM_LIT:50> ) <EOL> ml . write_input ( ) <EOL> m = ml . run_model ( silent = False ) <EOL> get_stp = [ <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ] <EOL> get_per = [ <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:1> , <NUM_LIT:1> , <NUM_LIT:1> ] <EOL> nswi_times = len ( get_per ) <EOL> zetafile = os . path . join ( dirs [ <NUM_LIT:0> ] , '<STR_LIT>' . format ( modelname ) ) <EOL> zobj = fu . CellBudgetFile ( zetafile ) <EOL> zeta = [ ] <EOL> for kk in zip ( get_stp , get_per ) : <EOL> zeta . append ( zobj . get_data ( kstpkper = kk , text = '<STR_LIT>' ) [ <NUM_LIT:0> ] ) <EOL> zeta = np . array ( zeta ) <EOL> nlay_swt = <NUM_LIT> <EOL> timprs = ( np . arange ( <NUM_LIT:8> ) + <NUM_LIT:1> ) * <NUM_LIT> <EOL> nprs = len ( timprs ) <EOL> ndecay = <NUM_LIT:4> <EOL> ibound = np . ones ( ( nlay_swt , nrow , ncol ) , '<STR_LIT:int>' ) <EOL> for k in range ( <NUM_LIT:0> , nlay_swt ) : <EOL> ibound [ k , <NUM_LIT:0> , ncol - <NUM_LIT:1> ] = - <NUM_LIT:1> <EOL> bot = np . zeros ( ( nlay_swt , nrow , ncol ) , np . float ) <EOL> zall = [ <NUM_LIT:0> , - <NUM_LIT> , - <NUM_LIT> , - <NUM_LIT> , - <NUM_LIT> , - <NUM_LIT> , - <NUM_LIT> ] <EOL> dz = <NUM_LIT> / nlay_swt <EOL> tb = np . arange ( nlay_swt ) * - dz - dz <EOL> sconc = np . zeros ( ( nlay_swt , nrow , ncol ) , np . float ) <EOL> icbund = np . ones ( ( nlay_swt , nrow , ncol ) , np . int ) <EOL> strt = np . zeros ( ( nlay_swt , nrow , ncol ) , np . float ) <EOL> pressure = <NUM_LIT:0.> <EOL> g = <NUM_LIT> <EOL> z = - dz / <NUM_LIT> <EOL> for k in range ( <NUM_LIT:0> , nlay_swt ) : <EOL> for i in range ( <NUM_LIT:0> , ncol ) : <EOL> bot [ k , <NUM_LIT:0> , i ] = tb [ k ] <EOL> if bot [ k , <NUM_LIT:0> , <NUM_LIT:0> ] >= - <NUM_LIT> : <EOL> sconc [ k , <NUM_LIT:0> , : ] = <NUM_LIT:0.> / <NUM_LIT> * <NUM_LIT> * <NUM_LIT> / <NUM_LIT> <EOL> else : <EOL> sconc [ k , <NUM_LIT:0> , : ] = <NUM_LIT> / <NUM_LIT> * <NUM_LIT> * <NUM_LIT> / <NUM_LIT> <EOL> icbund [ k , <NUM_LIT:0> , - <NUM_LIT:1> ] = - <NUM_LIT:1> <EOL> dense = <NUM_LIT> + <NUM_LIT> * sconc [ k , <NUM_LIT:0> , <NUM_LIT:0> ] <EOL> pressure += <NUM_LIT:0.5> * dz * dense * g <EOL> if k > <NUM_LIT:0> : <EOL> z = z - dz <EOL> denseup = <NUM_LIT> + <NUM_LIT> * sconc [ k - <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:0> ] <EOL> pressure += <NUM_LIT:0.5> * dz * denseup * g <EOL> strt [ k , <NUM_LIT:0> , : ] = z + pressure / dense / g <EOL> khb = ( <NUM_LIT> * <NUM_LIT> * <NUM_LIT> / <NUM_LIT> ) * <NUM_LIT> * <NUM_LIT> * <NUM_LIT> <EOL> kvb = ( <NUM_LIT> * <NUM_LIT> * <NUM_LIT> / <NUM_LIT> ) * <NUM_LIT> * <NUM_LIT> * <NUM_LIT> <EOL> ssb = <NUM_LIT> <EOL> sszb = <NUM_LIT> <EOL> kh = np . zeros ( ( nlay_swt , nrow , ncol ) , np . float ) <EOL> kv = np . zeros ( ( nlay_swt , nrow , ncol ) , np . float ) <EOL> ss = np . zeros ( ( nlay_swt , nrow , ncol ) , np . float ) <EOL> ssz = np . zeros ( ( nlay_swt , nrow , ncol ) , np . float ) <EOL> for k in range ( <NUM_LIT:0> , nlay_swt ) : <EOL> for i in range ( <NUM_LIT:0> , ncol ) : <EOL> f = r [ i ] * <NUM_LIT> * math . pi <EOL> kh [ k , <NUM_LIT:0> , i ] = khb * f <EOL> kv [ k , <NUM_LIT:0> , i ] = kvb * f <EOL> ss [ k , <NUM_LIT:0> , i ] = ssb * f <EOL> ssz [ k , <NUM_LIT:0> , i ] = sszb * f <EOL> itype = mt3 . Mt3dSsm . itype_dict ( ) <EOL> nwell = <NUM_LIT:1> <EOL> for k in range ( <NUM_LIT:0> , nlay_swt ) : <EOL> if bot [ k , <NUM_LIT:0> , <NUM_LIT:0> ] >= - <NUM_LIT> : <EOL> nwell = k + <NUM_LIT:1> <EOL> print ( '<STR_LIT>' . format ( nlay_swt , dz , nwell ) ) <EOL> well_data = { } <EOL> ssm_data = { } <EOL> wellQ = - <NUM_LIT> <EOL> wellbtm = - <NUM_LIT> <EOL> wellQpm = wellQ / abs ( wellbtm ) <EOL> for ip in range ( <NUM_LIT:0> , nper ) : <EOL> welllist = np . zeros ( ( nwell , <NUM_LIT:4> ) , np . float ) <EOL> ssmlist = [ ] <EOL> for iw in range ( <NUM_LIT:0> , nwell ) : <EOL> if ip == <NUM_LIT:0> : <EOL> q = wellQpm * dz <EOL> else : <EOL> q = <NUM_LIT:0.0> <EOL> welllist [ iw , <NUM_LIT:0> ] = iw <EOL> welllist [ iw , <NUM_LIT:1> ] = <NUM_LIT:0> <EOL> welllist [ iw , <NUM_LIT:2> ] = <NUM_LIT:0> <EOL> welllist [ iw , <NUM_LIT:3> ] = q <EOL> ssmlist . append ( [ iw , <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0.> , itype [ '<STR_LIT>' ] ] ) <EOL> well_data [ ip ] = welllist . copy ( ) <EOL> ssm_data [ ip ] = ssmlist <EOL> modelname = '<STR_LIT>' <EOL> swtexe_name = '<STR_LIT>' <EOL> if not skipRuns : <EOL> ml = mf . Modflow ( modelname , version = '<STR_LIT>' , exe_name = swtexe_name , model_ws = dirs [ <NUM_LIT:1> ] ) <EOL> discret = mf . ModflowDis ( ml , nrow = nrow , ncol = ncol , nlay = nlay_swt , delr = delr , delc = delc , <EOL> top = <NUM_LIT:0> , botm = bot , <EOL> laycbd = <NUM_LIT:0> , nper = nper , perlen = perlen , nstp = nstp , steady = True ) <EOL> bas = mf . ModflowBas ( ml , ibound = ibound , strt = strt ) <EOL> lpf = mf . ModflowLpf ( ml , hk = kh , vka = kv , ss = ss , sy = ssz , vkcb = <NUM_LIT:0> , laytyp = <NUM_LIT:0> , layavg = <NUM_LIT:1> ) <EOL> wel = mf . ModflowWel ( ml , stress_period_data = well_data ) <EOL> oc = mf . ModflowOc88 ( ml , save_head_every = <NUM_LIT> ) <EOL> pcg = mf . ModflowPcg ( ml , hclose = <NUM_LIT> , rclose = <NUM_LIT> , mxiter = <NUM_LIT:100> , iter1 = <NUM_LIT:50> ) <EOL> ml . write_input ( ) <EOL> mt = mt3 . Mt3dms ( modelname , '<STR_LIT>' , ml , model_ws = dirs [ <NUM_LIT:1> ] ) <EOL> adv = mt3 . Mt3dAdv ( mt , mixelm = - <NUM_LIT:1> , <EOL> percel = <NUM_LIT:0.5> , <EOL> nadvfd = <NUM_LIT:0> , <EOL> nplane = <NUM_LIT:4> , <EOL> mxpart = <NUM_LIT> , <EOL> itrack = <NUM_LIT:2> , <EOL> dceps = <NUM_LIT> , <EOL> npl = <NUM_LIT:16> , <EOL> nph = <NUM_LIT:16> , <EOL> npmin = <NUM_LIT:8> , <EOL> npmax = <NUM_LIT> ) <EOL> btn = mt3 . Mt3dBtn ( mt , icbund = icbund , prsity = ssz , ncomp = <NUM_LIT:1> , sconc = sconc , ifmtcn = - <NUM_LIT:1> , <EOL> chkmas = False , nprobs = <NUM_LIT:10> , nprmas = <NUM_LIT:10> , dt0 = <NUM_LIT:1.0> , ttsmult = <NUM_LIT:1.0> , <EOL> nprs = nprs , timprs = timprs , mxstrn = <NUM_LIT> ) <EOL> dsp = mt3 . Mt3dDsp ( mt , al = <NUM_LIT:0.> , trpt = <NUM_LIT:1.> , trpv = <NUM_LIT:1.> , dmcoef = <NUM_LIT:0.> ) <EOL> gcg = mt3 . Mt3dGcg ( mt , mxiter = <NUM_LIT:1> , iter1 = <NUM_LIT:50> , isolve = <NUM_LIT:1> , cclose = <NUM_LIT> ) <EOL> ssm = mt3 . Mt3dSsm ( mt , stress_period_data = ssm_data ) <EOL> mt . write_input ( ) <EOL> mswt = swt . Seawat ( modelname , '<STR_LIT>' , ml , mt , <EOL> exe_name = swtexe_name , model_ws = dirs [ <NUM_LIT:1> ] ) <EOL> vdf = swt . SeawatVdf ( mswt , iwtable = <NUM_LIT:0> , densemin = <NUM_LIT:0> , densemax = <NUM_LIT:0> , denseref = <NUM_LIT> , denseslp = <NUM_LIT> , firstdt = <NUM_LIT> ) <EOL> mswt . write_input ( ) <EOL> m = mswt . run_model ( silent = False ) <EOL> timprs <EOL> ucnfile = os . path . join ( dirs [ <NUM_LIT:1> ] , '<STR_LIT>' ) <EOL> uobj = fu . UcnFile ( ucnfile ) <EOL> times = uobj . get_times ( ) <EOL> print ( times ) <EOL> conc = np . zeros ( ( len ( times ) , nlay_swt , ncol ) , np . float ) <EOL> for idx , tt in enumerate ( times ) : <EOL> c = uobj . get_data ( totim = tt ) <EOL> for ilay in range ( <NUM_LIT:0> , nlay_swt ) : <EOL> for jcol in range ( <NUM_LIT:0> , ncol ) : <EOL> conc [ idx , ilay , jcol ] = c [ ilay , <NUM_LIT:0> , jcol ] <EOL> bot = np . zeros ( ( <NUM_LIT:1> , ncol , nlay ) , np . float ) <EOL> dz = <NUM_LIT> / float ( nlay - <NUM_LIT:1> ) <EOL> zall = - np . arange ( <NUM_LIT:0> , <NUM_LIT:100> + dz , dz ) <EOL> zall = np . append ( zall , - <NUM_LIT> ) <EOL> tb = - np . arange ( dz , <NUM_LIT:100> + dz , dz ) <EOL> tb = np . append ( tb , - <NUM_LIT> ) <EOL> for k in range ( <NUM_LIT:0> , nlay ) : <EOL> for i in range ( <NUM_LIT:0> , ncol ) : <EOL> bot [ <NUM_LIT:0> , i , k ] = tb [ k ] <EOL> swt_dz = <NUM_LIT> / nlay_swt <EOL> swt_tb = np . zeros ( ( nlay_swt ) , np . float ) <EOL> zc = - swt_dz / <NUM_LIT> <EOL> for klay in range ( <NUM_LIT:0> , nlay_swt ) : <EOL> swt_tb [ klay ] = zc <EOL> zc -= swt_dz <EOL> X , Z = np . meshgrid ( x , swt_tb ) <EOL> fwid , fhgt = <NUM_LIT> , <NUM_LIT> <EOL> flft , frgt , fbot , ftop = <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> <EOL> eps = <NUM_LIT> <EOL> lc = [ '<STR_LIT:r>' , '<STR_LIT:c>' , '<STR_LIT:g>' , '<STR_LIT:b>' , '<STR_LIT:k>' ] <EOL> cfig = [ '<STR_LIT:A>' , '<STR_LIT:B>' , '<STR_LIT:C>' , '<STR_LIT:D>' , '<STR_LIT:E>' , '<STR_LIT:F>' , '<STR_LIT>' , '<STR_LIT:H>' ] <EOL> inc = <NUM_LIT> <EOL> xsf = plt . figure ( figsize = ( fwid , fhgt ) , facecolor = '<STR_LIT:w>' ) <EOL> xsf . subplots_adjust ( wspace = <NUM_LIT> , hspace = <NUM_LIT> , left = flft , right = frgt , bottom = fbot , top = ftop ) <EOL> ax = xsf . add_subplot ( <NUM_LIT:4> , <NUM_LIT:2> , <NUM_LIT:1> ) <EOL> ax . text ( <NUM_LIT:0.0> , <NUM_LIT> , '<STR_LIT>' , transform = ax . transAxes , va = '<STR_LIT>' , ha = '<STR_LIT:left>' , size = '<STR_LIT>' ) <EOL> ax = xsf . add_subplot ( <NUM_LIT:4> , <NUM_LIT:2> , <NUM_LIT:2> ) <EOL> ax . text ( <NUM_LIT:0.0> , <NUM_LIT> , '<STR_LIT>' , transform = ax . transAxes , va = '<STR_LIT>' , ha = '<STR_LIT:left>' , size = '<STR_LIT>' ) <EOL> ax = xsf . add_subplot ( <NUM_LIT:4> , <NUM_LIT:2> , <NUM_LIT:1> ) <EOL> ax . plot ( [ - <NUM_LIT:1> , - <NUM_LIT:1> ] , [ - <NUM_LIT:1> , - <NUM_LIT:1> ] , '<STR_LIT>' , markersize = <NUM_LIT:3> , markeredgecolor = '<STR_LIT>' , markerfacecolor = '<STR_LIT:None>' , label = '<STR_LIT>' ) <EOL> ax . plot ( [ - <NUM_LIT:1> , - <NUM_LIT:1> ] , [ - <NUM_LIT:1> , - <NUM_LIT:1> ] , color = '<STR_LIT:k>' , linewidth = <NUM_LIT> , linestyle = '<STR_LIT>' , label = '<STR_LIT>' ) <EOL> ax . plot ( [ - <NUM_LIT:1> , - <NUM_LIT:1> ] , [ - <NUM_LIT:1> , - <NUM_LIT:1> ] , marker = '<STR_LIT:s>' , color = '<STR_LIT:k>' , linewidth = <NUM_LIT:0> , linestyle = '<STR_LIT:none>' , markeredgecolor = '<STR_LIT:w>' , <EOL> markerfacecolor = '<STR_LIT>' , label = '<STR_LIT>' ) <EOL> leg = ax . legend ( loc = '<STR_LIT>' , numpoints = <NUM_LIT:1> , ncol = <NUM_LIT:1> , labelspacing = <NUM_LIT:0.5> , borderaxespad = <NUM_LIT:1> , handlelength = <NUM_LIT:3> ) <EOL> leg . _drawFrame = False <EOL> for itime in range ( <NUM_LIT:0> , nswi_times ) : <EOL> zb = np . zeros ( ( ncol ) , np . float ) <EOL> zs = np . zeros ( ( ncol ) , np . float ) <EOL> for icol in range ( <NUM_LIT:0> , ncol ) : <EOL> for klay in range ( <NUM_LIT:0> , nlay ) : <EOL> ztop = float ( '<STR_LIT>' . format ( zall [ klay ] ) ) <EOL> zbot = float ( '<STR_LIT>' . format ( zall [ klay + <NUM_LIT:1> ] ) ) <EOL> zt = zeta [ itime , klay , <NUM_LIT:0> , icol ] <EOL> if ( ztop - zt ) > eps : <EOL> zs [ icol ] = zt <EOL> if itime < ndecay : <EOL> ic = itime <EOL> isp = ic * <NUM_LIT:2> + <NUM_LIT:1> <EOL> ax = xsf . add_subplot ( <NUM_LIT:4> , <NUM_LIT:2> , isp ) <EOL> else : <EOL> ic = itime - ndecay <EOL> isp = ( ic * <NUM_LIT:2> ) + <NUM_LIT:2> <EOL> ax = xsf . add_subplot ( <NUM_LIT:4> , <NUM_LIT:2> , isp ) <EOL> ax . text ( - <NUM_LIT> , <NUM_LIT> , cfig [ itime ] , transform = ax . transAxes , va = '<STR_LIT>' , ha = '<STR_LIT>' , size = '<STR_LIT>' ) <EOL> ax . plot ( x , zs , '<STR_LIT>' , markersize = <NUM_LIT:3> , markeredgecolor = '<STR_LIT>' , markerfacecolor = '<STR_LIT:None>' , label = '<STR_LIT>' ) <EOL> sc = ax . contour ( X , Z , conc [ itime , : , : ] , levels = [ <NUM_LIT> ] , colors = '<STR_LIT:k>' , linestyles = '<STR_LIT>' , linewidths = <NUM_LIT> , zorder = <NUM_LIT:30> ) <EOL> cc = ax . contourf ( X , Z , conc [ itime , : , : ] , levels = [ <NUM_LIT:0.0> , <NUM_LIT> , <NUM_LIT> ] , colors = [ '<STR_LIT:w>' , '<STR_LIT>' , '<STR_LIT:w>' ] ) <EOL> ax . set_xlim ( <NUM_LIT:0> , <NUM_LIT> ) <EOL> ax . set_ylim ( - <NUM_LIT:100> , - <NUM_LIT> ) <EOL> if itime < ndecay : <EOL> ax . set_ylabel ( '<STR_LIT>' ) <EOL> ax = xsf . add_subplot ( <NUM_LIT:4> , <NUM_LIT:2> , <NUM_LIT:7> ) <EOL> ax . set_xlabel ( '<STR_LIT>' ) <EOL> ax = xsf . add_subplot ( <NUM_LIT:4> , <NUM_LIT:2> , <NUM_LIT:8> ) <EOL> ax . set_xlabel ( '<STR_LIT>' ) <EOL> for itime in range ( <NUM_LIT:0> , nswi_times ) : <EOL> if itime < ndecay : <EOL> ic = itime <EOL> isp = ic * <NUM_LIT:2> + <NUM_LIT:1> <EOL> ax = xsf . add_subplot ( <NUM_LIT:4> , <NUM_LIT:2> , isp ) <EOL> else : <EOL> ic = itime - ndecay <EOL> isp = ( ic * <NUM_LIT:2> ) + <NUM_LIT:2> <EOL> ax = xsf . add_subplot ( <NUM_LIT:4> , <NUM_LIT:2> , isp ) <EOL> iyr = itime + <NUM_LIT:1> <EOL> if iyr > <NUM_LIT:1> : <EOL> ctxt = '<STR_LIT>' . format ( iyr ) <EOL> else : <EOL> ctxt = '<STR_LIT>' . format ( iyr ) <EOL> ax . text ( <NUM_LIT> , <NUM_LIT> , ctxt , transform = ax . transAxes , va = '<STR_LIT>' , ha = '<STR_LIT:right>' , size = '<STR_LIT>' ) <EOL> outfig = '<STR_LIT>' . format ( fext ) <EOL> xsf . savefig ( outfig , dpi = <NUM_LIT> ) <EOL> print ( '<STR_LIT>' , outfig ) </s>
