notebooks/PFLOTRAN_input_deck.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate PFLOTRAN input deck"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is used to generate the input deck that could be read by PFLOTRAN. It contains information such as numerical solver, grid discretization, time step and etc.\n",
    "\n",
    "Note: the `sink` and `cat` function is not directing output to file. The alternative is to redefine `cat` function ([post](https://stackoverflow.com/questions/51292852/how-to-use-sink-to-save-r-output-in-jupyter-notebook))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I/O files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-27T23:04:17.229422Z",
     "start_time": "2018-07-27T23:04:17.207Z"
    }
   },
   "outputs": [],
   "source": [
    "imodel=c(\"100x100x1\")\n",
    "model_dir = \"/global/project/projectdirs/m1800/pin/Reach_scale_model/\"\n",
    "## inputs\n",
    "fname_WellScreen = paste(model_dir, \"data/well_data/HEIS_300A_well_screen.csv\", sep = \"\")\n",
    "\n",
    "fname_material.h5 = c(\"HFR_material_river.h5\")\n",
    "fname_H.initial.h5 = c(\"HFR_H_Initial.h5\")\n",
    "# fname_H.BC.h5 = c(\"HFR_H_BC.h5\")\n",
    "\n",
    "fname_bc_dir = \"bc_6h_smooth/\"\n",
    "fname.DatumH = c(\"DatumH_Mass1_\")\n",
    "fname.Gradient = c(\"Gradients_Mass1_\")\n",
    "\n",
    "fname_mass_section = paste(model_dir, \"results/HFR_model_\", imodel, \"/mass_section.txt\", sep = \"\")\n",
    "## outputs\n",
    "fname_model_inputs.r = paste(model_dir, \"results/HFR_model_\", imodel, \"/model_inputs.r\", sep = \"\")\n",
    "fname_pflotran.in = paste(model_dir, \"Inputs/HFR_model_\", imodel, \"/pflotran_\", imodel, \"_new_iniH.in\", sep = \"\")\n",
    "\n",
    "is.flowBC = FALSE\n",
    "has.basalt = TRUE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-13T19:11:57.412148Z",
     "start_time": "2018-07-13T19:11:57.395Z"
    },
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "load(fname_model_inputs.r)\n",
    "# print(fname_pflotran.in)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## model parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-13T19:11:59.841987Z",
     "start_time": "2018-07-13T19:11:59.795Z"
    }
   },
   "outputs": [],
   "source": [
    "#-------------PARAMETERIZATION------------------#\n",
    "start.time = as.POSIXct(\"2007-03-28 12:00:00\",tz=\"GMT\",format=\"%Y-%m-%d %H:%M:%S\")\n",
    "end.time = as.POSIXct(\"2015-12-31 00:00:00\",tz=\"GMT\",format=\"%Y-%m-%d %H:%M:%S\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-13T19:12:01.265186Z",
     "start_time": "2018-07-13T19:12:01.251Z"
    }
   },
   "outputs": [],
   "source": [
    "# --------------------------times----------------------------------\n",
    "time.index = seq(from=start.time,to=end.time,by=\"1 hour\")\n",
    "ntime = length(time.index)\n",
    "nhours = ntime # from 2010-02-27 12:00:00 to 2017-07-01 00:00:00\n",
    "ini.ts = 0.01 #hour\n",
    "max.ts = 6 #hour"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-13T19:11:59.841987Z",
     "start_time": "2018-07-13T19:11:59.795Z"
    }
   },
   "outputs": [],
   "source": [
    "grid.n_pts = c(nx, ny, nz)\n",
    "grid.d_pts = c(idx,  idy,  idz)\n",
    "\n",
    "river_cond = 4.65e-13 # from mean(rand_2) conductance\n",
    "\n",
    "model_domain_origin = c(0,0,z0)\n",
    "\n",
    "a = 0 #rotation\n",
    "wall.stop = 47.8 # wallclock stop time, h\n",
    "check.pt = 720 # check point, h\n",
    "\n",
    "# solver\n",
    "ts.acceleration = 8\n",
    "max.ts.cuts = 20\n",
    "\n",
    "diff.coeff = 1e-9 # diffusion coefficient\n",
    "# recharge = 1.757e-9 #recharge rate 5.54 cm/yr from Fayer and Walters (1995)\n",
    "\n",
    "##----------------------- material -------------------------------##\n",
    "if (has.basalt) {\n",
    "    material.list = c(\"hanford\", \"cold_creek\", \"taylor_flats\", \"ringold_e\", \"ringold_lm\", \"ringold_a\", \"basalt\")\n",
    "} else {\n",
    "    material.list = c(\"hanford\", \"cold_creek\", \"taylor_flats\", \"ringold_e\", \"ringold_lm\", \"ringold_a\")\n",
    "}\n",
    "\n",
    "# material properties\n",
    "Khf = 7000 #m/d from 300A model, William's 2008\n",
    "hanford = list(ID = 1, porosity = 0.2, tortuosity = 1, perm.x = Khf/(9.5e11), perm.y = Khf/(9.5e11), perm.z = Khf/(9.5e12), cc = \"cc1\")\n",
    "# Kcc = 4.8 #m/d from Last et al\n",
    "Kcc = 100 #m/d from Paul's report\n",
    "cold_creek = list(ID = 2, porosity = 0.25, tortuosity = 1, perm.x = Kcc/(9.5e11), perm.y = Kcc/(9.5e11), perm.z = Kcc/(9.5e12), cc = \"cc1\")\n",
    "Ktf = 1 # m/d from franklin county report\n",
    "taylor_flats = list(ID = 3, porosity = 0.43, tortuosity = 1, perm.x = Ktf/(9.5e11), perm.y = Ktf/(9.5e11), perm.z = Ktf/(9.5e12), cc = \"cc2\")\n",
    "Kre = 40 # m/d from Paul's report\n",
    "ringold_e = list(ID = 4, porosity = 0.25, tortuosity = 1, perm.x = Kre/(9.5e11), perm.y = Kre/(9.5e11), perm.z = Kre/(9.5e12), cc = \"cc1\")\n",
    "Krlm = 1 #m/d from franklin\n",
    "ringold_lm = list(ID = 5, porosity = 0.43, tortuosity = 1, perm.x = Krlm/(9.5e11), perm.y = Krlm/(9.5e11), perm.z = Krlm/(9.5e12), cc = \"cc2\")\n",
    "Kra = 1 # similar to krlm\n",
    "ringold_a = list(ID = 6, porosity = 0.43, tortuosity = 1, perm.x = Kra/(9.5e11), perm.y = Kra/(9.5e11), perm.z = Kra/(9.5e12), cc = \"cc2\")\n",
    "\n",
    "Kb = 1e-3 # from Franklin\n",
    "basalt = list(ID = 9, porosity = 0.2, tortuosity = 1, perm.x = Kb/(9.5e11), perm.y = Kb/(9.5e11), perm.z = Kb/(9.5e12), cc = \"cc1\")\n",
    "\n",
    "K.list = list(hanford = hanford, cold_creek = cold_creek, taylor_flats = taylor_flats, ringold_e = ringold_e, ringold_lm = ringold_lm, \n",
    "              ringold_a = ringold_a, basalt =  basalt)\n",
    "# saturation functions\n",
    "cc.names = c(\"cc1\", \"cc2\")\n",
    "\n",
    "cc1 = list(M = 0.3391, alpha = 7.2727e-4, rsat = 0.16, max.cp = 1e8)\n",
    "cc2 = list(M = 0.7479, alpha = 1.4319e-4, rsat = 0.1299, max.cp = 1e8)\n",
    "\n",
    "cc.list = list(cc1 = cc1, cc2 = cc2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-13T19:12:00.916178Z",
     "start_time": "2018-07-13T19:12:00.906Z"
    }
   },
   "outputs": [],
   "source": [
    "# output opitons\n",
    "output.ts = 6 # hour\n",
    "h5.files = 73 # number of hdf5 files to write\n",
    "obs.ts = 6 # hour"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## add flow regions\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-11T00:20:06.515674Z",
     "start_time": "2018-07-11T00:20:06.462Z"
    }
   },
   "outputs": [],
   "source": [
    "##-----------------------flow region------------------------\n",
    "mass_section = read.delim(fname_mass_section, header = FALSE, stringsAsFactor = FALSE)\n",
    "slice.list = as.character(unlist(mass_section))\n",
    "\n",
    "# slice.list = as.character(seq(40, 332, 1)) \n",
    "river_region = paste(\"Mass1_\", slice.list, sep = \"\")\n",
    "flow_region = paste(\"Flow_Mass1_\", slice.list, sep = \"\")\n",
    "names(river_region)  = flow_region\n",
    "names(slice.list) = flow_region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-13T19:12:38.108440Z",
     "start_time": "2018-07-13T19:12:38.076Z"
    }
   },
   "outputs": [],
   "source": [
    "## ------------------- add solute tracer----------------------------\n",
    "solute.list = c(\"Solute_river\")\n",
    "\n",
    "## condition coupler, each flow region corresponds to one solute\n",
    "solute.coupler.list = c(rep(solute.list, length(flow_region)))\n",
    "names(solute.coupler.list)= flow_region\n",
    "\n",
    "## concentration list should have same number as of solute list\n",
    "concentration.list = c(\"Concentration_river\")\n",
    "names(concentration.list) = solute.list\n",
    "\n",
    "## tracer list should has same num. of solute list\n",
    "tracer.list = c(\"Tracer\")\n",
    "names(tracer.list) = concentration.list\n",
    "\n",
    "## conductance assign same num. as of flow region\n",
    "cond.list = c(rep(river_cond, length(flow_region)))\n",
    "\n",
    "names(cond.list) = flow_region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-13T19:12:39.327941Z",
     "start_time": "2018-07-13T19:12:39.277Z"
    }
   },
   "outputs": [],
   "source": [
    "## create river flow bc list\n",
    "\n",
    "fname.river_datum = paste(fname.DatumH, slice.list, \".txt\", sep = \"\")\n",
    "fname.river_gradient = paste(fname.Gradient, slice.list, \".txt\", sep = \"\")\n",
    "names(fname.river_datum)= flow_region\n",
    "names(fname.river_gradient)= flow_region\n",
    "\n",
    "# datasets\n",
    "flow.BC = c(\"Flow_West\", \"Flow_South\", \"Flow_North\", \"Flow_East\")\n",
    "# BC.all = c(\"BC_East\", \"BC_West\", \"BC_South\", \"BC_North\")\n",
    "BC.list = c(\"BC_West\", \"BC_South\", \"BC_North\", \"BC_East\")\n",
    "BC.face = c(\"West\", \"South\", \"North\", \"East\")\n",
    "names(flow.BC) = BC.face\n",
    "names(BC.list) = flow.BC\n",
    "max.buff = 200"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-13T19:12:40.770893Z",
     "start_time": "2018-07-13T19:12:40.761Z"
    }
   },
   "outputs": [],
   "source": [
    "# model corner coordinates\n",
    "xrange = c(0, xlen)\n",
    "yrange = c(0, ylen)\n",
    "zrange = c(z0, z0+zlen)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## obs wells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-13T19:12:42.763577Z",
     "start_time": "2018-07-13T19:12:42.698Z"
    }
   },
   "outputs": [],
   "source": [
    "#--------------- generate well obs input deck in Pflotran-----------------\n",
    "#read well screen info\n",
    "Wells = read.csv(fname_WellScreen, header = TRUE, stringsAsFactors=FALSE)\n",
    "Wells$WellName = paste(\"Well_\", Wells$WellName, sep=\"\")\n",
    "# Wells[,1] = paste(\"Well_\", Wells[,1], sep=\"\")\n",
    "rownames(Wells) = Wells$WellName\n",
    "# names(Wells) = c(\"wellname\", \"x\", \"y\", \"elev\", \"screen_top\", \"screen_bot\")\n",
    "Well = as.matrix(Wells[,3:7])\n",
    "\n",
    "\n",
    "Easting = Well[,1]\n",
    "Northing = Well[,2]\n",
    "elev = Well[,3]\n",
    "screen_top = Well[,4]\n",
    "screen_bot = Well[,5]\n",
    "\n",
    "\n",
    "# convert project coord to model coord, and rotate coordinates a-deg\n",
    "x = Easting -x0\n",
    "y = Northing -y0\n",
    "xx=x*cos(a)+y*sin(a)\n",
    "yy=y*cos(a)-x*sin(a)\n",
    "\n",
    "nwell = length(x)\n",
    "cells_z = seq((z0+0.5*grid.d_pts[3]), (z0+grid.n_pts[3]*grid.d_pts[3]), grid.d_pts[3])\n",
    "n_obs = nwell\n",
    "\n",
    "well.list = Wells$WellName\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# write input deck"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Redefine `cat` function so it directs output to external file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f <- file(fname_pflotran.in, open = \"wt\")\n",
    "cat <- function(...){\n",
    "    base::cat(..., file=f)\n",
    "    }\n",
    "# cat(\"Hello World\\n\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## beginning of file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-11T17:35:58.862747Z",
     "start_time": "2018-07-11T17:35:58.854Z"
    }
   },
   "outputs": [],
   "source": [
    "# sink(file = fname_pflotran.in, append = FALSE, type = \"output\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-13T19:13:16.327432Z",
     "start_time": "2018-07-13T19:13:16.191Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cat(\"## ===============BEGINNING of FILE=========================\")\n",
    "cat(\"\\n\")\n",
    "cat(\"SIMULATION\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  SIMULATION_TYPE SUBSURFACE\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  PROCESS_MODELS\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  SUBSURFACE_FLOW flow\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  MODE RICHARDS\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  /\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  SUBSURFACE_TRANSPORT transport\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  GLOBAL_IMPLICIT\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  /\")\n",
    "cat(\"\\n\")\n",
    "cat(\"/\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "cat(\"CHECKPOINT\")\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"  PERIODIC TIME\"), check.pt, \"h\")\n",
    "cat(\"\\n\")\n",
    "cat(\"/\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "cat(\"#  RESTART pflotran_bigplume-restart.chk\\n\")\n",
    "# use result from previoius simulaiton and reset time to zero\n",
    "cat(\"#  RESTART \\n\") \n",
    "cat(\"#  FILENAME pflotran_bigplume-restart.chk \\n\") \n",
    "cat(\"#  RESET_TO_TIME_ZERO\\n\")\n",
    "cat(\"##  REALIZATION_DEPENDENT\\n\")\n",
    "cat(\"#  END\\n\")\n",
    "\n",
    "cat(\"END\\n\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "cat(\"SUBSURFACE\")\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"WALLCLOCK_STOP\"), wall.stop, \"h\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## chemistry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** need to uncomment `Tracer_Age`, `DATABASE`, `ACTIVITY` if simulate age **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-12T21:11:46.828634Z",
     "start_time": "2018-07-12T21:11:46.780Z"
    }
   },
   "outputs": [],
   "source": [
    "cat('#=======================chemistry========================')\n",
    "cat(\"\\n\")\n",
    "cat(\"CHEMISTRY\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  PRIMARY_SPECIES\")\n",
    "cat(\"\\n\")\n",
    "for (itracer in tracer.list) {\n",
    "  cat(paste(\"   \", itracer))\n",
    "  cat(\"\\n\")\n",
    "}\n",
    "cat(\"#    Tracer_Age\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  /\")\n",
    "cat(\"\\n\")\n",
    "cat(\"#  DATABASE tracer_HFR.dat\")\n",
    "cat(\"\\n\")\n",
    "cat(\"#  ACTIVITY_COEFFICIENTS OFF\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  OUTPUT\")\n",
    "cat(\"\\n\")\n",
    "cat(\"    ALL\")\n",
    "cat(\"\\n\")\n",
    "cat(\"#    AGE\")\n",
    "cat(\"\\n\")\n",
    "cat(\"    TOTAL\")\n",
    "cat(\"\\n\")\n",
    "cat(\"/\")\n",
    "cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## solver options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-10T23:55:02.677283Z",
     "start_time": "2018-07-10T23:55:02.484Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cat('#=======================solver options========================')\n",
    "cat(\"\\n\")\n",
    "cat(\"TIMESTEPPER FLOW\")\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"  TS_ACCELERATION\", ts.acceleration))\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"  MAX_TS_CUTS\", max.ts.cuts))\n",
    "cat(\"\\n\")\n",
    "cat(\"#  MAX_STEPS 100\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "cat(\"\\n\")\n",
    "cat(\"TIMESTEPPER TRANSPORT\")\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"  TS_ACCELERATION\", ts.acceleration))\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"  MAX_TS_CUTS\", max.ts.cuts))\n",
    "cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "cat(\"NEWTON_SOLVER FLOW\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  VERBOSE_ERROR_MESSAGING\") ## add detailed error message\n",
    "cat(\"\\n\")\n",
    "cat(\"  MAXIT 20\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  RTOL 1.d-50\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  ATOL 1.d-50\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  STOL 1.d-60\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  ITOL_UPDATE 1.d0\")\n",
    "cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "cat(\"LINEAR_SOLVER FLOW\")\n",
    "cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "cat(\"NEWTON_SOLVER TRANSPORT\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  NO_INFINITY_NORM\")\n",
    "cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "cat(\"LINEAR_SOLVER TRANSPORT\")\n",
    "cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## discretization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`DXYZ` may not be used with `BOUNDS`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-13T20:09:00.797275Z",
     "start_time": "2018-07-13T20:09:00.642Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cat('#=======================discretization========================')\n",
    "cat(\"\\n\")\n",
    "cat(\"GRID\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  TYPE structured\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  NXYZ \")\n",
    "cat(grid.n_pts)\n",
    "cat(\"\\n\")\n",
    "cat(\"  ORIGIN \")\n",
    "cat(model_domain_origin)\n",
    "cat(\"\\n\")\n",
    "cat(\"  DXYZ \")\n",
    "\n",
    "cat(\"\\n\")\n",
    "cat(\"  \")\n",
    "cat(grid.d_pts[1])\n",
    "cat(\"\\n\")\n",
    "cat(\"  \")\n",
    "cat(grid.d_pts[2])\n",
    "cat(\"\\n\")\n",
    "cat(\"  \")\n",
    "cat(grid.d_pts[3])\n",
    "cat(\"\\n\")\n",
    "cat(\"/\")\n",
    "cat(\"\\n\")\n",
    "# cat(\"  BOUNDS\")\n",
    "# cat(\"\\n\")\n",
    "# cat(\"  \")\n",
    "# cat(xrange[1], yrange[1], zrange[1])\n",
    "# cat(\"\\n\")\n",
    "# cat(\"  \")\n",
    "# cat(xrange[2], yrange[2], zrange[2])\n",
    "# cat(\"\\n\")\n",
    "# cat(\"/\")\n",
    "# cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## fluid properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-10T23:57:01.974002Z",
     "start_time": "2018-07-10T23:57:01.928Z"
    }
   },
   "outputs": [],
   "source": [
    "cat('#=======================fluid properties========================')\n",
    "cat(\"\\n\")\n",
    "cat(\"FLUID_PROPERTY\")\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"  DIFFUSION_COEFFICIENT\", diff.coeff))\n",
    "cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-10T23:57:50.226278Z",
     "start_time": "2018-07-10T23:57:50.124Z"
    }
   },
   "outputs": [],
   "source": [
    "cat('#=======================datasets========================')\n",
    "cat(\"\\n\")\n",
    "\n",
    "if (is.flowBC) {\n",
    "    for (iBC in BC.list) {\n",
    "      cat(paste(\"DATASET\", iBC))\n",
    "      cat(\"\\n\")\n",
    "      cat(paste(\"  FILENAME\", fname_H.BC.h5))\n",
    "      cat(\"\\n\")\n",
    "      cat(paste(\"  HDF5_DATASET_NAME\", iBC))\n",
    "      cat(\"\\n\")\n",
    "      cat(paste(\"  MAX_BUFFER_SIZE\", max.buff))\n",
    "      cat(\"\\n\")\n",
    "      cat(\"END\")\n",
    "      cat(\"\\n\")\n",
    "      cat(\"\\n\")\n",
    "    }\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-10T23:57:50.226278Z",
     "start_time": "2018-07-10T23:57:50.124Z"
    }
   },
   "outputs": [],
   "source": [
    "cat(\"DATASET Initial_Head\")\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"  FILENAME\", fname_H.initial.h5))\n",
    "cat(\"\\n\")\n",
    "cat(\"  HDF5_DATASET_NAME Initial_Head\")\n",
    "cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "cat(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## material properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-10T23:58:24.162906Z",
     "start_time": "2018-07-10T23:58:24.099Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cat('#=======================material properties========================')\n",
    "cat(\"\\n\")\n",
    "for (imaterial in material.list) {\n",
    "  \n",
    "  cat(paste(\"MATERIAL_PROPERTY\", imaterial))\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"  ID\"), K.list[[imaterial]]$ID)\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"  POROSITY\"), K.list[[imaterial]]$porosity)\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"  TORTUOSITY\"), K.list[[imaterial]]$tortuosity)\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"  CHARACTERISTIC_CURVES\", K.list[[imaterial]]$cc))\n",
    "  cat(\"\\n\")\n",
    "  cat(\"  PERMEABILITY\")\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"  PERM_X\"), K.list[[imaterial]]$perm.x)\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"  PERM_Y\"), K.list[[imaterial]]$perm.y)\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"  PERM_Z\"), K.list[[imaterial]]$perm.z)\n",
    "  cat(\"\\n\")\n",
    "  cat('  /')\n",
    "  cat(\"\\n\")\n",
    "  cat('END')\n",
    "  cat(\"\\n\")\n",
    "  cat(\"\\n\")\n",
    "  \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## saturation functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-10T23:59:05.544355Z",
     "start_time": "2018-07-10T23:59:05.456Z"
    }
   },
   "outputs": [],
   "source": [
    "cat('#=======================saturation fuctions========================')\n",
    "cat(\"\\n\")\n",
    "for (icc in cc.names) {\n",
    "  \n",
    "  cat(paste(\"CHARACTERISTIC_CURVES\", icc))\n",
    "  cat(\"\\n\")\n",
    "  cat(\"  SATURATION_FUNCTION VAN_GENUCHTEN\")\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"    M\"), cc.list[[icc]]$M)\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"    ALPHA\"), cc.list[[icc]]$alpha)\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"    LIQUID_RESIDUAL_SATURATION\"), cc.list[[icc]]$rsat)\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"    MAX_CAPILLARY_PRESSURE\"), cc.list[[icc]]$max.cp)\n",
    "  cat(\"\\n\")\n",
    "  cat('  /')\n",
    "  cat(\"\\n\")\n",
    "  cat(\"  PERMEABILITY_FUNCTION BURDINE_VG_LIQ\")\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"    M\"), cc.list[[icc]]$M)\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"    LIQUID_RESIDUAL_SATURATION\"), cc.list[[icc]]$rsat)\n",
    "  cat(\"\\n\")\n",
    "  cat('  /')\n",
    "  cat(\"\\n\")\n",
    "  cat('END')\n",
    "  cat(\"\\n\")\n",
    "  cat(\"\\n\")\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## output options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-17T23:33:54.938727Z",
     "start_time": "2018-07-17T23:33:54.885Z"
    }
   },
   "outputs": [],
   "source": [
    "cat(\"\\n\")\n",
    "cat('#=======================output options========================')\n",
    "cat(\"\\n\")\n",
    "cat(\"OUTPUT\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  VARIABLES\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  LIQUID_PRESSURE\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  LIQUID_HEAD\\n\")\n",
    "cat(\"  LIQUID_SATURATION\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  MATERIAL_ID_KLUDGE_FOR_VISIT\")  # add composite material ids\n",
    "cat(\"\\n\")\n",
    "cat(\"  RESIDUAL\\n\")\n",
    "cat('  /')\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  SNAPSHOT_FILE\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  NO_PRINT_INITIAL\")\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"#  PERIODIC TIME\"), output.ts, \"h\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  PERIODIC TIME 6 h between 32400 h and 76800 h\")\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"#  FORMAT HDF5 MULTIPLE_FILES TIMES_PER_FILE\"), h5.files)\n",
    "cat(\"\\n\")\n",
    "cat(\"  FORMAT HDF5 SINGLE_FILE\")\n",
    "cat(\"\\n\")\n",
    "cat('  /')\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  MASS_BALANCE_FILE\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  NO_PRINT_INITIAL\")\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"#  PERIODIC TIME\"), output.ts, \"h\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  PERIODIC TIME 6 h between 32400 h and 76800 h\")\n",
    "cat(\"\\n\")\n",
    "cat('  /')\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  OBSERVATION_FILE\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  NO_PRINT_INITIAL\\n\")\n",
    "cat(\"  PERIODIC TIME 6 h between 32400 h and 76800 h\")\n",
    "cat(\"\\n\")\n",
    "cat('  /')\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  VELOCITY_AT_CENTER \") # OUTPUT darcy velocity at cell center\n",
    "cat(\"\\n\")\n",
    "cat(\"  VELOCITY_AT_FACE \")  # output darcy velocity at cell face\n",
    "cat(\"\\n\")\n",
    "cat('END')\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat('#=======================times========================')\n",
    "cat(\"\\n\")\n",
    "cat(\"TIME\")\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"  FINAL_TIME\", nhours, \"h\"))\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"  INITIAL_TIMESTEP_SIZE\", ini.ts, \"h\"))\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"  MAXIMUM_TIMESTEP_SIZE\", max.ts, \"h\"))\n",
    "cat(\"\\n\")\n",
    "cat('END')\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-11T00:03:09.398812Z",
     "start_time": "2018-07-11T00:03:09.164Z"
    }
   },
   "outputs": [],
   "source": [
    "cat('#=======================regions========================')\n",
    "cat(\"\\n\")\n",
    "cat('REGION all')\n",
    "cat('\\n')\n",
    "cat('  COORDINATES')\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[1], yrange[1], zrange[1])\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[2], yrange[2], zrange[2])\n",
    "cat('\\n')\n",
    "cat('  /')\n",
    "cat('\\n')\n",
    "cat('END')\n",
    "cat('\\n')\n",
    "cat('\\n')\n",
    "\n",
    "cat('REGION Bottom')\n",
    "cat('\\n')\n",
    "cat('FACE BOTTOM')\n",
    "cat('\\n')\n",
    "cat('  COORDINATES')\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[1], yrange[1], zrange[1])\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[2], yrange[2], zrange[1])\n",
    "cat('\\n')\n",
    "cat('  /')\n",
    "cat('\\n')\n",
    "cat('END')\n",
    "cat('\\n')\n",
    "cat('\\n')\n",
    "\n",
    "cat('REGION Top')\n",
    "cat('\\n')\n",
    "cat('FACE TOP')\n",
    "cat('\\n')\n",
    "cat('  COORDINATES')\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[1], yrange[1], zrange[2])\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[2], yrange[2], zrange[2])\n",
    "cat('\\n')\n",
    "cat('  /')\n",
    "cat('\\n')\n",
    "cat('END')\n",
    "cat('\\n')\n",
    "cat('\\n')\n",
    "\n",
    "cat('REGION South')\n",
    "cat('\\n')\n",
    "cat('FACE SOUTH')\n",
    "cat('\\n')\n",
    "cat('  COORDINATES')\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[1], yrange[1], zrange[1])\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[2], yrange[1], zrange[2])\n",
    "cat('\\n')\n",
    "cat('  /')\n",
    "cat('\\n')\n",
    "cat('END')\n",
    "cat('\\n')\n",
    "cat('\\n')\n",
    "\n",
    "cat('REGION North')\n",
    "cat('\\n')\n",
    "cat('FACE NORTH')\n",
    "cat('\\n')\n",
    "cat('  COORDINATES')\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[1], yrange[2], zrange[1])\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[2], yrange[2], zrange[2])\n",
    "cat('\\n')\n",
    "cat('  /')\n",
    "cat('\\n')\n",
    "cat('END')\n",
    "cat('\\n')\n",
    "cat('\\n')\n",
    "\n",
    "cat('REGION West')\n",
    "cat('\\n')\n",
    "cat('FACE WEST')\n",
    "cat('\\n')\n",
    "cat('  COORDINATES')\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[1], yrange[1], zrange[1])\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[1], yrange[2], zrange[2])\n",
    "cat('\\n')\n",
    "cat('  /')\n",
    "cat('\\n')\n",
    "cat('END')\n",
    "cat('\\n')\n",
    "cat('\\n')\n",
    "\n",
    "cat('REGION East')\n",
    "cat('\\n')\n",
    "cat('FACE EAST')\n",
    "cat('\\n')\n",
    "cat('  COORDINATES')\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[2], yrange[1], zrange[1])\n",
    "cat('\\n')\n",
    "cat(\"  \")\n",
    "cat(xrange[2], yrange[2], zrange[2])\n",
    "cat('\\n')\n",
    "cat('  /')\n",
    "cat('\\n')\n",
    "cat('END')\n",
    "cat('\\n')\n",
    "cat('\\n')\n",
    "\n",
    "# cat('REGION Inj_Well')\n",
    "# cat('\\n')\n",
    "# cat('  COORDINATES')\n",
    "# cat('\\n')\n",
    "# cat(xx[inj_wellname], yy[inj_wellname], screen_bot[inj_wellname])\n",
    "# cat('\\n')\n",
    "# cat(xx[inj_wellname], yy[inj_wellname], screen_top[inj_wellname])\n",
    "# cat('\\n')\n",
    "# cat('  /')\n",
    "# cat('\\n')\n",
    "# cat('/')\n",
    "# cat('\\n')\n",
    "\n",
    "# %%%%%%%%%%%%%%%%%%%river regions%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
    "\n",
    "for (iriver in river_region) {\n",
    "  cat(paste('REGION ', iriver, sep = '' ))\n",
    "  cat('\\n')\n",
    "  cat(paste('  FILE ', fname_material.h5))\n",
    "  cat('\\n')\n",
    "  cat('END')\n",
    "  cat(\"\\n\")\n",
    "  cat(\"\\n\")\n",
    "}\n",
    "cat(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### add obs wells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-11T18:48:43.087557Z",
     "start_time": "2018-07-11T18:48:42.883Z"
    }
   },
   "outputs": [],
   "source": [
    "#------------- %put observation wells in \"regions\" of input deck----------------------\n",
    "cat('# Observation wells')\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "for (iwell in well.list) {\n",
    "  \n",
    "  if (screen_top[iwell] >= z0 & screen_bot[iwell] <= z0+zlen & (screen_top[iwell] - screen_bot[iwell]) >= dz[1]) {\n",
    "    z = cells_z[which((cells_z >= screen_bot[iwell]) & (cells_z <= screen_top[iwell]))]\n",
    "    nz = length(z)\n",
    "    # wellname = rownames(Well)[iwell]\n",
    "    for (j in 1:nz) {\n",
    "      cat(paste('REGION ', iwell, '_', j, sep = ''))\n",
    "      cat('\\n')\n",
    "      cat('  COORDINATE ')\n",
    "      cat(xx[iwell],yy[iwell], z[j])\n",
    "      cat('\\n')\n",
    "      cat('/')\n",
    "      cat('\\n')\n",
    "      cat('\\n')\n",
    "    }\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## flow conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-11T00:06:25.833741Z",
     "start_time": "2018-07-11T00:06:25.728Z"
    }
   },
   "outputs": [],
   "source": [
    "cat(\"#==================flow conditions=========================\")\n",
    "cat(\"\\n\")\n",
    "cat(\"MINIMUM_HYDROSTATIC_PRESSURE -1.d0\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "cat(\"FLOW_CONDITION Initial\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  TYPE\")\n",
    "cat(\"\\n\")\n",
    "cat(\"    PRESSURE hydrostatic\")\n",
    "cat(\"\\n\")\n",
    "cat('  /')\n",
    "cat(\"\\n\")\n",
    "cat(\"  CYCLIC\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  DATUM DATASET Initial_Head\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  PRESSURE 101325\")\n",
    "cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "# cat(\"FLOW_CONDITION Recharge\")\n",
    "# cat(\"\\n\")\n",
    "# cat(\"  TYPE\")\n",
    "# cat(\"\\n\")\n",
    "# cat(\"    FLUX neumann\")\n",
    "# cat(\"\\n\")\n",
    "# cat('/')\n",
    "# cat(\"\\n\")\n",
    "# cat(paste(\"    FLUX\"), recharge)\n",
    "# cat(\"\\n\")\n",
    "# cat('/')\n",
    "# cat(\"\\n\")\n",
    "\n",
    "\n",
    "#%%%%%%%%%%%%%%%%%%%% river regions %%%%%%%%%%%%%%%%%%%%%%%%%\n",
    "for (iflow in flow_region) {\n",
    "  cat(paste(\"FLOW_CONDITION\"), iflow)\n",
    "  cat(\"\\n\")\n",
    "  cat(\"  INTERPOLATION LINEAR\\n\")\n",
    "  cat(\"  TYPE\")\n",
    "  cat(\"\\n\")\n",
    "  cat(\"  PRESSURE conductance\")\n",
    "  cat(\"\\n\")\n",
    "  cat('  /')\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"  CONDUCTANCE\"), cond.list[iflow])\n",
    "  cat(\"\\n\")\n",
    "  cat(\"  CYCLIC\")\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"  DATUM file \", fname_bc_dir, fname.river_datum[iflow], sep = \"\"))\n",
    "  cat(\"\\n\")\n",
    "  cat(\"  PRESSURE 101325\")\n",
    "  cat(\"\\n\")\n",
    "  cat(\"  GRADIENT\")\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"  PRESSURE file \",fname_bc_dir, fname.river_gradient[iflow], sep = \"\"))\n",
    "  cat(\"\\n\")\n",
    "  cat('  /')\n",
    "  cat(\"\\n\")\n",
    "  cat('/')\n",
    "  cat(\"\\n\")\n",
    "  cat(\"\\n\")\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-11T00:06:25.833741Z",
     "start_time": "2018-07-11T00:06:25.728Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## flow BCs\n",
    "\n",
    "if (is.flowBC) {\n",
    "    for (iflow.BC in flow.BC) {\n",
    "      cat(paste(\"FLOW_CONDITION\"), iflow.BC)\n",
    "      cat(\"\\n\")\n",
    "      cat(\"  TYPE\")\n",
    "      cat(\"\\n\")\n",
    "      cat(\"  PRESSURE hydrostatic\")\n",
    "      cat(\"\\n\")\n",
    "      cat(\"  /\")\n",
    "      cat(\"\\n\")\n",
    "      cat(\"  CYCLIC\")\n",
    "      cat(\"\\n\")\n",
    "      cat(paste(\"  DATUM DATASET\"), BC.list[iflow.BC])\n",
    "      cat(\"\\n\")\n",
    "      cat(\"  PRESSURE 101325\")\n",
    "      cat(\"\\n\")\n",
    "      cat(\"/\")\n",
    "      cat(\"\\n\")\n",
    "      cat(\"\\n\")\n",
    "    }\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## transport conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-11T00:06:55.392341Z",
     "start_time": "2018-07-11T00:06:55.292Z"
    }
   },
   "outputs": [],
   "source": [
    "cat(\"#==================transport conditions=========================\")\n",
    "cat(\"\\n\")\n",
    "cat(\"TRANSPORT_CONDITION Initial\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  TYPE dirichlet_zero_gradient\") \n",
    "cat(\"\\n\")\n",
    "cat(\"  CONSTRAINT_LIST\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  \")\n",
    "cat(paste(0, \"Concentration_initial\"))\n",
    "cat(\"\\n\")\n",
    "cat(\"/\")\n",
    "cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "for (isolute in solute.list) {\n",
    "  cat(paste(\"TRANSPORT_CONDITION\"), isolute)\n",
    "  cat(\"\\n\")\n",
    "  cat(\"  TYPE dirichlet_zero_gradient\") # use dirichlet zero gradient when mass is moving in and out of boundary\n",
    "  cat(\"\\n\")\n",
    "  cat(\"  CONSTRAINT_LIST\")\n",
    "  cat(\"\\n\")\n",
    "  cat(\"  \")\n",
    "  cat(paste(0, concentration.list[isolute]))\n",
    "  cat(\"\\n\")\n",
    "  cat(\"/\")\n",
    "  cat(\"\\n\")\n",
    "  cat(\"END\")\n",
    "  cat(\"\\n\")\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** need to comment out `Tracer_Age` for age simulation **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-11T00:07:27.455589Z",
     "start_time": "2018-07-11T00:07:27.392Z"
    }
   },
   "outputs": [],
   "source": [
    "cat(\"#==================constraints=========================\")\n",
    "cat(\"\\n\")\n",
    "for (iconc in concentration.list) {\n",
    "  cat(paste(\"CONSTRAINT\", iconc))\n",
    "  cat(\"\\n\")\n",
    "  cat(\"  CONCENTRATIONS\")\n",
    "  cat(\"\\n\")\n",
    "  for (itracer in tracer.list) {\n",
    "    if (itracer == tracer.list[iconc]) {\n",
    "      cat(\"   \")\n",
    "      cat(paste(itracer), \"1.d0 T\")\n",
    "    } else {\n",
    "      cat(\"   \")\n",
    "      cat(paste(itracer), \"1.d-20 T\")\n",
    "    }\n",
    "    cat(\"\\n\")\n",
    "\n",
    "  }\n",
    "  cat(\"#    Tracer_Age  1.d-20 T\")\n",
    "  cat(\"\\n\")  \n",
    "  cat(\"  /\")\n",
    "  cat(\"\\n\")\n",
    "  cat(\"END\")\n",
    "  cat(\"\\n\")\n",
    "  cat(\"\\n\")\n",
    "}\n",
    "\n",
    "cat(\"\\n\")\n",
    "cat(\"CONSTRAINT Concentration_initial\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  CONCENTRATIONS\")\n",
    "cat(\"\\n\")\n",
    "for (itracer in tracer.list) {\n",
    "  cat(\"   \")\n",
    "  cat(paste(itracer), \"1.d-20 T\")\n",
    "  cat(\"\\n\")\n",
    "}\n",
    "cat(\"#    Tracer_Age  1.d-20 T\")\n",
    "cat(\"\\n\")  \n",
    "cat(\"  /\")\n",
    "cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "cat(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## condition couplers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-11T00:08:35.225518Z",
     "start_time": "2018-07-11T00:08:35.148Z"
    }
   },
   "outputs": [],
   "source": [
    "cat(\"#==================condition couplers=========================\")\n",
    "# \n",
    "cat(\"\\n\")\n",
    "cat(\"INITIAL_CONDITION\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  FLOW_CONDITION Initial\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  TRANSPORT_CONDITION Initial\")\n",
    "cat(\"\\n\")\n",
    "cat(\"  REGION all\")\n",
    "cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")\n",
    "\n",
    "#%%%%%%%%%%%%%%%%%%%%% River BC coupling %%%%%%%%%%%%%%%%%%%%%%%%%%\n",
    "for (iflow in flow_region) {\n",
    "  cat(paste(\"BOUNDARY_CONDITION River_\", slice.list[iflow], sep = \"\"))\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"  FLOW_CONDITION\"), iflow)\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"  TRANSPORT_CONDITION\"), solute.coupler.list[iflow])\n",
    "  cat(\"\\n\")\n",
    "  cat(paste(\"  REGION\"), river_region[iflow])\n",
    "  cat(\"\\n\")\n",
    "  cat(\"END\")\n",
    "  cat(\"\\n\")\n",
    "  cat(\"\\n\")\n",
    "}\n",
    "\n",
    "# for (iface in BC.face) {\n",
    "#   cat(paste(\"BOUNDARY_CONDITION\", iface))\n",
    "#   cat(\"\\n\")\n",
    "#   cat(paste(\"  FLOW_CONDITION\", flow.BC[iface]))\n",
    "#   cat(\"\\n\")\n",
    "#   if (iface == \"North\") {\n",
    "#     cat(paste(\"  TRANSPORT_CONDITION\", solute.list[5]))\n",
    "#   } else {\n",
    "#     cat(\"  TRANSPORT_CONDITION Initial\")\n",
    "#   }\n",
    "#   cat(\"\\n\")\n",
    "#   cat(paste(\"  REGION\"), iface) \n",
    "#   cat(\"\\n\")\n",
    "#   cat(\"END\")\n",
    "#   cat(\"\\n\")\n",
    "#   cat(\"\\n\")\n",
    "# }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-11T00:08:35.225518Z",
     "start_time": "2018-07-11T00:08:35.148Z"
    }
   },
   "outputs": [],
   "source": [
    "if (is.flowBC) {\n",
    "    \n",
    "    for (iface in BC.face) {\n",
    "      cat(paste(\"BOUNDARY_CONDITION\", iface))\n",
    "      cat(\"\\n\")\n",
    "      cat(paste(\"  FLOW_CONDITION\", flow.BC[iface]))\n",
    "      cat(\"\\n\")\n",
    "        cat(\"  TRANSPORT_CONDITION Initial\")\n",
    "      cat(\"\\n\")\n",
    "      cat(paste(\"  REGION\"), iface)\n",
    "      cat(\"\\n\")\n",
    "      cat(\"END\")\n",
    "      cat(\"\\n\")\n",
    "      cat(\"\\n\")\n",
    "    }\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## stratigraphy couplers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-27T23:04:23.259086Z",
     "start_time": "2018-07-27T23:04:23.238Z"
    }
   },
   "outputs": [],
   "source": [
    "cat(\"#==================stratigraphy couplers=========================\")\n",
    "cat(\"\\n\")\n",
    "cat(\"STRATA\")\n",
    "cat(\"\\n\")\n",
    "cat(paste(\"  FILE\", fname_material.h5))\n",
    "cat(\"\\n\")\n",
    "cat(\"END\")\n",
    "cat(\"\\n\")\n",
    "cat(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## observation points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-11T00:09:49.482730Z",
     "start_time": "2018-07-11T00:09:49.235Z"
    }
   },
   "outputs": [],
   "source": [
    "cat(\"#==================observation points=========================\")\n",
    "cat('\\n')\n",
    "# % putting in \"observation points\" in input deck\n",
    "for (iwell in well.list){\n",
    "  \n",
    "  if (screen_top[iwell] >= z0 & screen_bot[iwell] <= z0+zlen & (screen_top[iwell] - screen_bot[iwell]) >= dz[1]) {\n",
    "    z = cells_z[which(cells_z >= screen_bot[iwell] & (cells_z <= screen_top[iwell]))]\n",
    "    nz = length(z)\n",
    "    # wellname = rownames(Well)[iwell]\n",
    "    for (j in 1:nz){\n",
    "      \n",
    "      cat('OBSERVATION')\n",
    "      cat('\\n')\n",
    "      cat(paste('  REGION ', iwell, '_', j, sep = ''))\n",
    "      cat('\\n')\n",
    "      cat('  VELOCITY ')\n",
    "      cat('\\n')\n",
    "      cat('/')\n",
    "      cat('\\n')\n",
    "      cat('\\n')\n",
    "    }\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## end of file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-11T00:10:06.181871Z",
     "start_time": "2018-07-11T00:10:06.169Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cat(\"#==================END=========================\")\n",
    "cat(\"\\n\")\n",
    "cat(\"END_SUBSURFACE\")\n",
    "# sink()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-11T17:35:58.862747Z",
     "start_time": "2018-07-11T17:35:58.854Z"
    }
   },
   "outputs": [],
   "source": [
    "close(f)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.5.1"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "oldHeight": 650,
   "position": {
    "height": "40px",
    "left": "11px",
    "right": "20px",
    "top": "606px",
    "width": "367px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "varInspector_section_display": "none",
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}