## 
! {RMODFLOW} is still in its experimental lifecycle stage.
## 
! Use at your own risk, and submit issues here:
## 
! <https://github.com/rogiersbart/RMODFLOW/issues>

RMODFLOW provides rmf_create_*, rmf_read_* and rmf_write_* functions for individual MODFLOW packages. In most cases however, it is useful to read, write and create entire MODFLOW models. For this, RMODFLOW provides the top-level rmf_create, rmf_read and rmf_write functions. The main advantage of using these top-level functions, besides the obvious less verbose code, is that all options are set correctly: free/fixed formats, binary, additional arguments etc. are all set correctly by these functions.

Read

To read in an entire model, use rmf_read() on the MODFLOW name file. This will read in all the input packages active during the simulation. The only additional argument that has to be specified by the user is the precision of binary files (“single” by default).

rma <- rmf_example_file("rocky-mountain-arsenal.nam") %>%
  rmf_read()
#> ---------------------------------
#>   Reading NAM package from file C:/Users/brogiers/AppData/Local/Temp/Rtmp8kwVIi/temp_libpath30a852a11ff6/RMODFLOW/extdata/rocky-mountain-arsenal.nam
#> ---------------------------------
#>   Reading DIS package from file C:/Users/brogiers/AppData/Local/Temp/Rtmp8kwVIi/temp_libpath30a852a11ff6/RMODFLOW/extdata/rocky-mountain-arsenal.dis
#> ---------------------------------
#>   Reading BAS6 package from file C:/Users/brogiers/AppData/Local/Temp/Rtmp8kwVIi/temp_libpath30a852a11ff6/RMODFLOW/extdata/rocky-mountain-arsenal.bas
#> ---------------------------------
#>   Reading LPF package from file C:/Users/brogiers/AppData/Local/Temp/Rtmp8kwVIi/temp_libpath30a852a11ff6/RMODFLOW/extdata/rocky-mountain-arsenal.lpf
#> ---------------------------------
#>   Reading WEL package from file C:/Users/brogiers/AppData/Local/Temp/Rtmp8kwVIi/temp_libpath30a852a11ff6/RMODFLOW/extdata/rocky-mountain-arsenal.wel
#> ---------------------------------
#>   Reading CHD package from file C:/Users/brogiers/AppData/Local/Temp/Rtmp8kwVIi/temp_libpath30a852a11ff6/RMODFLOW/extdata/rocky-mountain-arsenal.chd
#> ---------------------------------
#>   Reading PCG package from file C:/Users/brogiers/AppData/Local/Temp/Rtmp8kwVIi/temp_libpath30a852a11ff6/RMODFLOW/extdata/rocky-mountain-arsenal.pcg
#> ---------------------------------
#>   Reading OC package from file C:/Users/brogiers/AppData/Local/Temp/Rtmp8kwVIi/temp_libpath30a852a11ff6/RMODFLOW/extdata/rocky-mountain-arsenal.oc
#> ---------------------------------
#>   Reading LMT6 package from file C:/Users/brogiers/AppData/Local/Temp/Rtmp8kwVIi/temp_libpath30a852a11ff6/RMODFLOW/extdata/rocky-mountain-arsenal.lmt

To also include output files, set output = TRUE but be warned: this might become slow if output files are large. To suppress the writing of information to the console, set verbose = FALSE.

rmf_read returns a list of class modflow which holds all the packages. It can be useful to “unpack” this list with list2env(model, .Globalenv).

Create

To create a RMODFLOW modflow object, use rmf_create(). This function takes as input RMODFLOW objects of class rmf_package or a single list with rmf_package objects. If a nam object is not provided, it is created (silently). To overwrite the cell-by-cell flow package number of all objects that have a this number, specify a cbc value.

modflow <- rmf_create(rma$dis, rma$bas, rma$lpf, rma$oc, rma$pcg, rma$wel, rma$chd, cbc = 88)
str(modflow)
#> List of 8
#>  $ dis:List of 16
#>   ..$ nlay  : num 1
#>   ..$ nrow  : num 54
#>   ..$ ncol  : num 46
#>   ..$ nper  : num 1
#>   ..$ itmuni: num 1
#>   ..$ lenuni: num 2
#>   ..$ laycbd: num 0
#>   ..$ delr  : num [1:46] 99.8 99.8 99.8 99.8 99.8 ...
#>   ..$ delc  : num [1:54] 100 100 100 100 100 ...
#>   ..$ top   : 'rmf_2d_array' num [1:54, 1:46] 70.7 69.3 68 66.6 65.2 ...
#>   .. ..- attr(*, "dimlabels")= chr [1:2] "i" "j"
#>   ..$ botm  : 'rmf_3d_array' num [1:54, 1:46, 1] 58.7 57.3 56 54.6 53.2 ...
#>   .. ..- attr(*, "dimlabels")= chr [1:3] "i" "j" "k"
#>   ..$ perlen: num 6.31e+08
#>   ..$ nstp  : num 1
#>   ..$ tsmult: num 1
#>   ..$ sstr  : chr "SS"
#>   ..$ prj   :List of 3
#>   .. ..$ origin  : Named num [1:3] 1696 3049 0
#>   .. .. ..- attr(*, "names")= chr [1:3] "x" "y" "z"
#>   .. ..$ rotation: num 0
#>   .. ..$ crs     :List of 2
#>   .. .. ..$ input: chr NA
#>   .. .. ..$ wkt  : chr NA
#>   .. .. ..- attr(*, "class")= chr "crs"
#>   .. ..- attr(*, "class")= chr "prj"
#>   ..- attr(*, "comment")= chr [1:6] " Discretization File created on 2017-05-30 by ModelMuse version 3.9.0.0." " Upper left corner: (1696.20102745308, 8473.29309713687)" " Lower left corner: (1696.20102745308, 3049.15053472316)" " Upper right corner: (6289.17627630449, 8473.29309713687)" ...
#>   ..- attr(*, "class")= chr [1:2] "dis" "rmf_package"
#>  $ bas:List of 10
#>   ..$ xsection    : logi FALSE
#>   ..$ chtoch      : logi TRUE
#>   ..$ free        : logi TRUE
#>   ..$ printtime   : logi TRUE
#>   ..$ showprogress: logi FALSE
#>   ..$ stoperror   : logi FALSE
#>   ..$ stoper      : num NA
#>   ..$ ibound      : 'rmf_3d_array' int [1:54, 1:46, 1] 0 0 0 0 0 0 0 0 0 0 ...
#>   .. ..- attr(*, "dimlabels")= chr [1:3] "i" "j" "k"
#>   ..$ hnoflo      : num -1e+20
#>   ..$ strt        : 'rmf_3d_array' num [1:54, 1:46, 1] 0 0 0 0 0 0 0 0 0 0 ...
#>   .. ..- attr(*, "dimlabels")= chr [1:3] "i" "j" "k"
#>   ..- attr(*, "comment")= chr " Basic Package file created on 2017-05-30 by ModelMuse version 3.9.0.0."
#>   ..- attr(*, "class")= chr [1:2] "bas" "rmf_package"
#>  $ lpf:List of 17
#>   ..$ ilpfcb            : num 88
#>   ..$ hdry              : num -2e+20
#>   ..$ nplpf             : num 0
#>   ..$ storagecoefficient: logi FALSE
#>   ..$ constantcv        : logi FALSE
#>   ..$ thickstrt         : logi FALSE
#>   ..$ nocvcorrection    : logi FALSE
#>   ..$ novfc             : logi FALSE
#>   ..$ noparcheck        : logi FALSE
#>   ..$ laytyp            : num 0
#>   ..$ layavg            : num 0
#>   ..$ chani             : num -1
#>   ..$ layvka            : num 0
#>   ..$ laywet            : num 0
#>   ..$ hk                : 'rmf_3d_array' num [1:54, 1:46, 1] 1e-04 1e-04 1e-04 1e-04 1e-04 1e-04 1e-04 1e-04 1e-04 1e-04 ...
#>   .. ..- attr(*, "dimlabels")= chr [1:3] "i" "j" "k"
#>   ..$ hani              : 'rmf_3d_array' num [1:54, 1:46, 1] 1 1 1 1 1 1 1 1 1 1 ...
#>   .. ..- attr(*, "dimlabels")= chr [1:3] "i" "j" "k"
#>   ..$ vka               : 'rmf_3d_array' num [1:54, 1:46, 1] 1e-05 1e-05 1e-05 1e-05 1e-05 1e-05 1e-05 1e-05 1e-05 1e-05 ...
#>   .. ..- attr(*, "dimlabels")= chr [1:3] "i" "j" "k"
#>   ..- attr(*, "comment")= chr " LPF: Layer Property Flow package file created on 2017-05-30 by ModelMuse version 3.9.0.0."
#>   ..- attr(*, "class")= chr [1:2] "lpf" "rmf_package"
#>  $ oc :List of 22
#>   ..$ ibouun        : logi NA
#>   ..$ cboufm        : logi NA
#>   ..$ iddnun        : num 38
#>   ..$ cddnfm        : chr "(10(1X1PE13.5))"
#>   ..$ iddnfm        : logi NA
#>   ..$ ihedun        : num 37
#>   ..$ chedfm        : chr "(10(1X1PE13.5))"
#>   ..$ ihedfm        : logi NA
#>   ..$ ibound_label  : logi FALSE
#>   ..$ drawdown_label: logi TRUE
#>   ..$ head_label    : logi TRUE
#>   ..$ aux           : logi TRUE
#>   ..$ compact_budget: logi TRUE
#>   ..$ print_head    : logi FALSE
#>   ..$ print_drawdown: logi FALSE
#>   ..$ save_head     : logi TRUE
#>   ..$ save_drawdown : logi TRUE
#>   ..$ save_ibound   : logi FALSE
#>   ..$ iperoc        : num 1
#>   ..$ itsoc         : num 1
#>   ..$ save_budget   : logi TRUE
#>   ..$ print_budget  : logi TRUE
#>   ..- attr(*, "comment")= chr " Output Control file created on 2017-05-30 by ModelMuse version 3.9.0.0."
#>   ..- attr(*, "class")= chr [1:2] "oc" "rmf_package"
#>  $ pcg:List of 11
#>   ..$ mxiter  : num 20
#>   ..$ iter1   : num 30
#>   ..$ npcond  : num 1
#>   ..$ ihcofadd: num 0
#>   ..$ hclose  : num 0.001
#>   ..$ rclose  : num 1000
#>   ..$ relax   : num 1
#>   ..$ nbpol   : num 1
#>   ..$ iprpcg  : num 1
#>   ..$ mutpcg  : num 0
#>   ..$ damppcg : num 1
#>   ..- attr(*, "comment")= chr " PCG: Preconditioned Conjugate Gradient package file created on 2017-05-30 by ModelMuse version 3.9.0.0."
#>   ..- attr(*, "class")= chr [1:2] "pcg" "rmf_package"
#>  $ wel:List of 10
#>   ..$ np       : num 0
#>   ..$ mxl      : num 0
#>   ..$ instances: NULL
#>   ..$ mxact    : num 4
#>   ..$ itmp     : Named num 4
#>   .. ..- attr(*, "names")= chr "list_1"
#>   ..$ iwelcb   : num 88
#>   ..$ option   : Named logi FALSE
#>   .. ..- attr(*, "names")= chr "noprint"
#>   ..$ aux      : chr "IFACE"
#>   ..$ data     :Classes 'rmf_list' and 'data.frame': 4 obs. of  7 variables:
#>   .. ..$ k        : int [1:4] 1 1 1 1
#>   .. ..$ i        : int [1:4] 13 13 42 44
#>   .. ..$ j        : int [1:4] 29 30 15 23
#>   .. ..$ q        : num [1:4] 0.0125 0.0125 -0.001 -0.002
#>   .. ..$ IFACE    : num [1:4] 0 0 0 0
#>   .. ..$ parameter: logi [1:4] FALSE FALSE FALSE FALSE
#>   .. ..$ name     : chr [1:4] "list_1" "list_1" "list_1" "list_1"
#>   ..$ kper     :'data.frame':    1 obs. of  2 variables:
#>   .. ..$ kper  : int 1
#>   .. ..$ list_1: logi TRUE
#>   ..- attr(*, "class")= chr [1:2] "wel" "rmf_package"
#>   ..- attr(*, "comment")= chr " WEL: Well package file created on 2017-05-30 by ModelMuse version 3.9.0.0."
#>  $ chd:List of 9
#>   ..$ np       : num 0
#>   ..$ mxl      : num 0
#>   ..$ instances: NULL
#>   ..$ mxact    : num 122
#>   ..$ itmp     : Named num 122
#>   .. ..- attr(*, "names")= chr "list_1"
#>   ..$ option   : Named logi FALSE
#>   .. ..- attr(*, "names")= chr "noprint"
#>   ..$ aux      : chr "IFACE"
#>   ..$ data     :Classes 'rmf_list' and 'data.frame': 122 obs. of  8 variables:
#>   .. ..$ k        : int [1:122] 1 1 1 1 1 1 1 1 1 1 ...
#>   .. ..$ i        : int [1:122] 1 2 1 2 1 2 3 1 2 3 ...
#>   .. ..$ j        : int [1:122] 8 8 9 9 10 10 10 11 11 11 ...
#>   .. ..$ shead    : num [1:122] 75 75 75 75 75 75 75 75 75 75 ...
#>   .. ..$ ehead    : num [1:122] 75 75 75 75 75 75 75 75 75 75 ...
#>   .. ..$ IFACE    : num [1:122] 0 0 0 0 0 0 0 0 0 0 ...
#>   .. ..$ parameter: logi [1:122] FALSE FALSE FALSE FALSE FALSE FALSE ...
#>   .. ..$ name     : chr [1:122] "list_1" "list_1" "list_1" "list_1" ...
#>   ..$ kper     :'data.frame':    1 obs. of  2 variables:
#>   .. ..$ kper  : int 1
#>   .. ..$ list_1: logi TRUE
#>   ..- attr(*, "class")= chr [1:2] "chd" "rmf_package"
#>   ..- attr(*, "comment")= chr " CHD: Time-Variant Specified-Head package file created on 2017-05-30 by ModelMuse version 3.9.0.0."
#>  $ nam:Classes 'nam', 'rmf_package' and 'data.frame':    11 obs. of  4 variables:
#>   ..$ ftype  : chr [1:11] "LIST" "DIS" "BAS6" "LPF" ...
#>   ..$ nunit  : num [1:11] 700 701 702 703 704 705 706 707 37 38 ...
#>   ..$ fname  : chr [1:11] "output.lst" "input.dis" "input.bas" "input.lpf" ...
#>   ..$ options: chr [1:11] "REPLACE" NA NA NA ...
#>   ..- attr(*, "dir")= chr ""
#>  - attr(*, "class")= chr "modflow"

All MODFLOW simulations must have at least a DIS, BAS, solver and flow package. If any of those are missing or duplicated, rmf_create will throw an error.

modflow <- rmf_create(rma$dis, rma$bas, rma$lpf, rma$oc)
#> Error: Please specify at least a dis, bas, solver and flow package. Only specify one of each.

Write

To write a model, use rmf_write() and pass a modflow object and the path to a MODFLOW name file.

rmf_write(modflow, file = "input.nam", verbose = FALSE)

You can exclude certain packages from writing by supplying an exclude argument with their RMODFLOW abbreviations. This can be useful if you e.g. want try out excluding a package from the simulation without recreating the entire modflow object. Because this bypasses the creation of the modflow object, you can exclude required packages such as the solver or the DIS file without RMODFLOW erroring out, so be careful.

rmf_write(modflow, file = "input.nam", verbose = FALSE, exclude = c("wel", "chd"))