diff --git a/Manifest.toml b/Manifest.toml
index 65fb7b8e..4ebf660c 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -39,12 +39,36 @@ git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
 uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
 version = "0.5.10"
 
+[[Blosc]]
+deps = ["Blosc_jll"]
+git-tree-sha1 = "84cf7d0f8fd46ca6f1b3e0305b4b4a37afe50fd6"
+uuid = "a74b3585-a348-5f62-a45c-50e91977d574"
+version = "0.7.0"
+
+[[Blosc_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Lz4_jll", "Pkg", "Zlib_jll", "Zstd_jll"]
+git-tree-sha1 = "e747dac84f39c62aff6956651ec359686490134e"
+uuid = "0b7ba130-8d10-5ba8-a3d6-c5182647fed9"
+version = "1.21.0+0"
+
+[[BufferedStreams]]
+deps = ["Compat", "Test"]
+git-tree-sha1 = "5d55b9486590fdda5905c275bb21ce1f0754020f"
+uuid = "e1450e63-4bb3-523b-b2a4-4ffa8c0fd77d"
+version = "1.0.0"
+
 [[Bzip2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
 uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
 version = "1.0.8+0"
 
+[[CFTime]]
+deps = ["Dates", "Printf"]
+git-tree-sha1 = "bca6cb6ee746e6485ca4535f6cc29cf3579a0f20"
+uuid = "179af706-886a-5703-950a-314cd64e0468"
+version = "0.1.1"
+
 [[CSV]]
 deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
 git-tree-sha1 = "c0a735698d1a0a388c5c7ae9c7fb3da72fd5424e"
@@ -299,9 +323,9 @@ uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
 
 [[GLFW_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"]
-git-tree-sha1 = "dba1e8614e98949abfa60480b13653813d8f0157"
+git-tree-sha1 = "0c603255764a1fa0b61752d2bec14cfbd18f7fe8"
 uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
-version = "3.3.5+0"
+version = "3.3.5+1"
 
 [[GR]]
 deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
@@ -356,6 +380,12 @@ git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
 uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
 version = "1.0.2"
 
+[[HDF5]]
+deps = ["Blosc", "Compat", "HDF5_jll", "Libdl", "Mmap", "Random", "Requires"]
+git-tree-sha1 = "83173193dc242ce4b037f0263a7cc45afb5a0b85"
+uuid = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+version = "0.15.6"
+
 [[HDF5_jll]]
 deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "Libdl", "OpenSSL_jll", "Pkg", "Zlib_jll"]
 git-tree-sha1 = "fd83fa0bde42e01952757f01149dd968c06c4dba"
@@ -567,6 +597,18 @@ git-tree-sha1 = "491a883c4fef1103077a7f648961adbf9c8dd933"
 uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
 version = "2.1.2"
 
+[[Lz4_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "5d494bc6e85c4c9b626ee0cab05daa4085486ab1"
+uuid = "5ced341a-0733-55b8-9ab6-a4889d929147"
+version = "1.9.3+0"
+
+[[MAT]]
+deps = ["BufferedStreams", "CodecZlib", "HDF5", "SparseArrays"]
+git-tree-sha1 = "5c62992f3d46b8dce69bdd234279bb5a369db7d5"
+uuid = "23992714-dd62-5051-b70f-ba57cb901cac"
+version = "0.10.1"
+
 [[MacroTools]]
 deps = ["Markdown", "Random"]
 git-tree-sha1 = "5a5bc6bf062f0f95e62d0fe0a2d99699fed82dd9"
@@ -604,6 +646,12 @@ uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 [[MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
 
+[[NCDatasets]]
+deps = ["CFTime", "DataStructures", "Dates", "NetCDF_jll", "Printf"]
+git-tree-sha1 = "5da406d9624f25909a6f556bd8d5c1deaa189ee6"
+uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
+version = "0.11.7"
+
 [[NLSolversBase]]
 deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"]
 git-tree-sha1 = "144bab5b1443545bc4e791536c9f1eacb4eed06a"
diff --git a/Project.toml b/Project.toml
index 6759618f..a10c32b8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,6 +17,8 @@ Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
+NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
diff --git a/README.md b/README.md
index 17a08ea1..d1c70bdd 100644
--- a/README.md
+++ b/README.md
@@ -6,17 +6,14 @@
 [![Coverage](https://codecov.io/gh/ggebbie/TMI.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/ggebbie/TMI.jl)
 
 * Total Matrix Intercomparison codes for Julia\
-G Jake Gebbie, WHOI, ggebbie@whoi.edu 
+Started by G Jake Gebbie, WHOI, ggebbie@whoi.edu 
 
-# References 
-
-Gebbie, G., and P. Huybers:  "Total matrix intercomparison: A method for resolving the geometry of water mass pathways", J. Phys. Oceanogr., 40(8), doi:10.1175/2010JPO4272.1, 1710-1728, 2010. 
-
-Gebbie, G., and P. Huybers. "How is the ocean filled?", Geophys. Res. Lett., 38, L06604, doi:10.1029/2011GL046769, 2011 
+* See the function list in the documentation linked through the badge above
 
-Gebbie, G., and P. Huybers, "The mean age of ocean waters inferred from radiocarbon observations", 2012, JPO.
+* The MATLAB version of the code is in maintenance mode and is available at https://github.com/ggebbie/TMI 
 
-Gebbie, G., "How much did Glacial North Atlantic Water shoal?", 2014, Paleoceanography.
+* After setting up the environment (instructions below), check that all tests pass via the following shell command in the repository base directory:
+`julia --project=@. test/runtests.jl`
 
 # Requirements
 
@@ -34,7 +31,7 @@ Restart the REPL.\
 `Conda.add("shapely",channel="conda-forge")`\
 `Conda.add("cartopy",channel="conda-forge")`\
 
-This should set up cartopy v. 0.20.0 which can download coastlines from Amazon Web Services.
+This should set up cartopy v0.20.0 which can download coastlines from Amazon Web Services.
 
 # Setting up project environment
 
@@ -71,15 +68,23 @@ An example:\
 
 # Data files
 
-The Julia code is designed to download input files from Google Drive and to place them in the `data` directory. If that doesn't work, extract data from Google Drive using your favorite method, such as the script `readTMIfromGoogleDrive.sh` at a bash shell prompt. 
+The Julia code is designed to download input files from Google Drive and to place them in the `data` directory. If that doesn't work, extract data from Google Drive using your favorite method or download manually at: https://drive.google.com/drive/folders/1nFDSINbst84pK68aWwRGBVfYZkfN1mUR?usp=sharing . Available TMI versions include:
 
-`/bin/sh readTMIfromGoogleDrive.sh`
+`modern_90x45x33_GH10_GH12` : TMI version with 4x4 degree horizontal
+                  resolution and 33 levels  (G & H 2010), \
+				  Includes the input data from the WGHC (Gouretski & Koltermann 2005) 
+ 
+`modern_180x90x33_GH10_GH12` : TMI version with 2x2 degree horizontal
+                  resolution and 33 levels  (G & H 2010), \
+				  Includes the input data from the WGHC (Gouretski & Koltermann 2005) 
 
-Or download manually at: https://drive.google.com/file/d/1Zycnx6_nifRrJo8XWMdlCFv4ODBpi-i7/view?usp=sharing .
+`modern_90x45x33_unpub12` : TMI version with 4x4 degree horizontal
+                  resolution and 33 levels  (unpublished 2012), \
+				  Includes a steady-state climatology of global tracers
 
-TMI_4deg_2010.mat : TMI version with 4x4 degree horizontal
-                  resolution and 33 levels  (G & H 2010), \
-				  Includes TMI climatology of ocean properties 
+`modern_90x45x33_G14` : TMI version with 4x4 degree horizontal
+                  resolution and 33 levels  (Gebbie 2014), \
+				  Doesn't rely upon a bottom spreading parameterization and solves for mixed-layer depth
 
 # Functions
 
@@ -97,9 +102,9 @@ Steps:
 1. Use PkgTemplates to make git repo.\
 	2. new empty repository on GitHub.\
     3. Then push an existing repository from the command line:
-    `git remote add origin git@github.com:ggebbie/TMI.jl.git`
-    `git branch -M main`
-    `git push -u origin main`
+    `git remote add origin git@github.com:ggebbie/TMI.jl.git`\
+    `git branch -M main`\
+    `git push -u origin main`\
 
 4. Run the following Julia code
 
@@ -118,9 +123,9 @@ Steps:
         Documenter{GitHubActions}(),
         Develop(),
     ],
-             )`
+             )`\
 
-`t("TMI.jl")`
+`t("TMI.jl")`\
 
 # MATLAB version of code
 
@@ -141,3 +146,13 @@ Version 6.2, July 2015, added sq.m function,
                         fixed d_all to properly divide Atlantic/Pacific and put White Sea into Arctic.\
 Version 7, Sept. 2016, major improvements to transient run: 2 types of initial conditions and boundary conditions.\
 Version 8, Jan. 2021, bug fixes, especially those found by Elaine McDonagh
+
+# References 
+
+Gebbie, G., and P. Huybers:  "Total matrix intercomparison: A method for resolving the geometry of water mass pathways", J. Phys. Oceanogr., 40(8), doi:10.1175/2010JPO4272.1, 1710-1728, 2010. 
+
+Gebbie, G., and P. Huybers. "How is the ocean filled?", Geophys. Res. Lett., 38, L06604, doi:10.1029/2011GL046769, 2011 
+
+Gebbie, G., and P. Huybers, "The mean age of ocean waters inferred from radiocarbon observations", 2012, JPO.
+
+Gebbie, G., "How much did Glacial North Atlantic Water shoal?", 2014, Paleoceanography.
diff --git a/scripts/convert_tmi_to_julia.jl b/scripts/convert_tmi_to_julia.jl
deleted file mode 100644
index 6c1e3fa1..00000000
--- a/scripts/convert_tmi_to_julia.jl
+++ /dev/null
@@ -1,232 +0,0 @@
-using MAT
-using LinearAlgebra
-using SparseArrays
-
-# load circulation file
-Afile = "/home/gebbie/Dropbox/mcode/TMI_v8/A_4deg_2010.mat"
-Avars = matread(Afile)
-
-tracerfile = "/home/gebbie/Dropbox/mcode/TMI_v8/tracerobs_4deg_33lev_TMI.mat"
-tracervars = matread(tracerfile)
-
-it = convert(Array{Int,1},vec(Avars["i"]));
-jt = convert(Array{Int,1},vec(Avars["j"]));
-kt = convert(Array{Int,1},vec(Avars["k"]));
-Itmi = CartesianIndex.(it,jt,kt)
-
-function initfld(vector::Array{Float64,2},I::Array{CartesianIndex{3},1})
-# function for changing TMI vector to field
-
-    #convert2vec!(vector)
-    vector = convert(Array{Float64,1},vec(vector))
-    println(size(vector))
-    nv = length(vector)
-
-    nx = Tuple(I[1])[1]
-    ny = Tuple(I[1])[2]
-    nz = Tuple(I[1])[3]
-    
-    for nn = 2:nv
-        nx = max(nx,Tuple(I[nn])[1])
-        ny = max(ny,Tuple(I[nn])[2])
-        nz = max(nz,Tuple(I[nn])[3])
-    end
-
-    field = NaN.*zeros(nx,ny,nz)
-    #- a comprehension
-    [field[I[n]]=vector[n] for n ∈ 1:nv]
-
-    return field
-end
-
-function convert2vec!(vector::Array{Float64,2})
-    vector = convert(Array{Float64,1},vec(vector))
-    return vector
-end
-
-function convert2intvec!(vector::Array{Float64,2})
-    convert(Array{Int64,1},vec(vector))
-end
-function convert2int!(vector::Array{Float64,2})
-    convert(Array{Int64,2},vector)
-end
-
-# update names and types in dictionary
-TMIgrids = Dict("lon" => convert2intvec!(Avars["LON"]),
-                "lat" => convert2intvec!(Avars["LAT"]),
-                "depth" => convert2intvec!(Avars["DEPTH"]))
-#                "Itmi" => Itmi,
-#                 "A" => Avars["A"])
-                 
-TMIfields = Dict("qpo4" => initfld(Avars["dP"],Itmi),
-                "θ" => initfld(tracervars["Tobs"],Itmi),
-                "σθ" => initfld(tracervars["Terr"],Itmi),
-                "sp" => initfld(tracervars["Sobs"],Itmi),
-                "σsp" => initfld(tracervars["Serr"],Itmi),
-                "δ18ow" => initfld(tracervars["O18obs"],Itmi),
-                "σδ18ow" => initfld(tracervars["O18err"],Itmi),
-                "po4" => initfld(tracervars["Pobs"],Itmi),
-                "σpo4" => initfld(tracervars["Perr"],Itmi),
-                "no3" => initfld(tracervars["Nobs"],Itmi),
-                "σno3" => initfld(tracervars["Nerr"],Itmi),
-                "o2" => initfld(tracervars["Oobs"],Itmi),
-                "σo2" => initfld(tracervars["Oerr"],Itmi),
-                "δ13c" => initfld(tracervars["C13obs"],Itmi),
-                "σδ13c" => initfld(tracervars["C13err"],Itmi))
-
-## write to NetCDF
-filenetcdf = "TMI_4deg_2010.nc"
-
-TMIgridsatts = Dict("lon" => Dict("longname" => "Longitude", "units" => "degrees east"),
-                    "lat" => Dict("longname" => "Latitude", "units" => "degrees north"),
-                    "depth" => Dict("longname" => "depth", "units" => "m"))
-
-TMIfieldsatts = Dict("θ" => Dict("longname" => "potential temperature", "units" => "°C"),
-                     "σθ" => Dict("longname" => "1σ standard error in potential temperature", "units" => "°C"),
-                     "sp" => Dict("longname" => "practical salinity", "units" => "PSS-78"),
-                     "σsp" => Dict("longname" => "1σ standard error in practical salinity", "units" => "PSS-78"),
-                     "δ18ow" => Dict("longname" => "oxygen-18 to oxygen-16 ratio in seawater", "units" => "‰ VSMOW"),
-                     "σδ18ow" => Dict("longname" => "1σ standard error in oxygen-18 to oxygen-16 ratio in seawater", "units" => "‰ VSMOW"),
-                     "po4" => Dict("longname" => "phosphate", "units" => "μmol/kg"),
-                     "σpo4" => Dict("longname" => "1σ standard error in phosphate", "units" => "μmol/kg"),
-                     "qpo4" => Dict("longname" => "local source of phosphate", "units" => "μmol/kg"),
-                     "no3" => Dict("longname" => "nitrate", "units" => "μmol/kg"),
-                     "σno3" => Dict("longname" => "1σ standard error in nitrate", "units" => "μmol/kg"),
-                     "o2" => Dict("longname" => "dissolved oxygen", "units" => "μmol/kg"),
-                     "σo2" => Dict("longname" => "1σ standard error in dissolved oxygen", "units" => "μmol/kg"),
-                     "δ13c" => Dict("longname" => "carbon-13 to carbon-12 ratio in DIC", "units" => "‰ PDB"),
-                     "σδ13c" => Dict("longname" => "1σ standard error in carbon-13 to carbon-12 ratio in DIC", "units" => "‰ PDB"))
-            
-
-# make new netcdf file.
-isfile(filenetcdf) && rm(filenetcdf)
-
-# iterate in TMIgrids Dictionary to write to NetCDF.
-for (varname,varvals) in TMIfields
-  println(varname)
-  nccreate(filenetcdf,varname,"lon",lon,TMIgridsatts["lon"],"lat",lat,TMIgridsatts["lat"],"depth",depth,TMIgridsatts["depth"],atts=TMIfieldsatts[varname])
-  ncwrite(varvals, filenetcdf,varname)
-end
-
-isrc = A_CSC[2]
-idst = A_CSC[1]
-ival = A_CSC[3]
-
-# add the circulation matrix: problem can't store sparse matrix.
-varname= "m"
-faceatts = Dict("longname" => "TMI grid face number")
-matts =  Dict("longname" => "TMI water-mass (sparse) matrix values","version"=> "2010")
-nccreate(filenetcdf,varname,"face",1:nface,faceatts,atts=matts)
-ncwrite(ival, filenetcdf,varname)
-
-varname= "i"
-destatts = Dict("longname" => "gridcell number of destination (row value)")
-nccreate(filenetcdf,varname,"face",1:nface,faceatts,atts=destatts)
-ncwrite(idst, filenetcdf,varname)
-
-varname= "j"
-sourceatts = Dict("longname" => "gridcell number of source (column value)")
-nccreate(filenetcdf,varname,"face",1:nface,faceatts,atts=sourceatts)
-ncwrite(isrc, filenetcdf,varname)
-
-varname = "xgrid"
-xgridatts = Dict("longname" => "gridcell Cartesian x-location")
-traceratts = Dict("longname" => "gridcell index")
-nccreate(filenetcdf,varname,"index number",1:nv,traceratts,atts=xgridatts)
-ncwrite(it, filenetcdf,varname)
-
-varname = "ygrid"
-ygridatts = Dict("longname" => "gridcell Cartesian y-location")
-traceratts = Dict("longname" => "gridcell index")
-nccreate(filenetcdf,varname,"index number",1:nv,traceratts,atts=ygridatts)
-ncwrite(jt, filenetcdf,varname)
-
-varname = "zgrid"
-zgridatts = Dict("longname" => "gridcell Cartesian z-location")
-traceratts = Dict("longname" => "gridcell index")
-nccreate(filenetcdf,varname,"index number",1:nv,traceratts,atts=zgridatts)
-ncwrite(kt, filenetcdf,varname)
-
-
-##############################################################
-## end of 4-jan-2021 work
-##############################################################
-
-
-function extract(d)
-           expr = quote end
-           for (k, v) in d
-               push!(expr.args, :($(Symbol(k)) = $v))
-           end
-           eval(expr)
-           return
-       end
-
-extract(Dict_TMI)
-
- dry = isnan.(po4)
- wet = .!dry
-
-# get po4fld
- @time po4fld = NaN.*zeros(nx,ny,nz);
- @time tmivec2fld!(po4fld,po4,indextmi);
-
-# get po4vector: slower process.
- @time po4vec = NaN.*zeros(nv);
- @time tmifld2vec!(po4vec,po4fld,indextmi);
-
-
-function tmivec2fld!(field::Array{Real,3},vector::Array{Real,1},index::Array{CartesianIndex{3},1})
-
-    nv = length(index)
-    #- a comprehension
-    [field[index[n]]=vector[n] for n ∈ 1:nv]
-
-    return field
-end
-
-function tmifld2vec!(vector::Array{Real,1},field::Array{Real,3},index::Array{CartesianIndex{3},1})
-
-    nv = length(index)
-    
-    #- a comprehension
-    [vector[n] = field[indextmi[n]] for n ∈ 1:nv];
-    return vector
-end
-
-
-Alu = lu(A)
-nfield = size(A,1)
-
-
-
-# vector to field via a comprehension (slow: 0.1 s)
-# @time po4fld_test = NaN.*zeros(nx,ny,nz);
-# @time [po4fld_test[indextmi[n]]=po4[n] for n ∈ 1:nv];
-
-# field to vector via a comprehension
-# @time po4 = [po4fld[indextmi[n]] for n ∈ 1:nv];
-#- also works -# @time [po4[n] = po4fld[indextmi[n]] for n ∈ 1:nv];
-
-#file = matopen(Afile)
-#for (key, values) in Avars
-#    println(key); #print(value)
-    #stringer = $key=value
-#    eval($key=value)
-#end
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/scripts/ex1.trackpathways.jl b/scripts/ex1.trackpathways.jl
index c72115f9..c8b75f54 100644
--- a/scripts/ex1.trackpathways.jl
+++ b/scripts/ex1.trackpathways.jl
@@ -13,7 +13,8 @@ using TMI, BenchmarkTools, PyPlot, PyCall
 
 pygui(true) #needed for Atom, not sure what it will do in other places
 
-TMIversion = "TMI_2010_2012_4x4x33"
+TMIversion = "modern_90x45x33_GH10_GH12"
+A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)
 
 #- define the surface patch by the bounding latitude and longitude.
 latbox = [50,60]; # 50 N -> 60 N, for example.
@@ -22,7 +23,7 @@ latbox = [50,60]; # 50 N -> 60 N, for example.
 lonbox = [-50,0]; # 50 W -> prime meridian
 
 # do numerical analysis
-c,γ = trackpathways(TMIversion,latbox,lonbox)
+c = trackpathways(Alu,latbox,lonbox,γ)
 
 # do plotting (could be a function)
 plotextent(latbox, lonbox)
diff --git a/scripts/ex2.watermassdistribution.jl b/scripts/ex2.watermassdistribution.jl
new file mode 100644
index 00000000..33c193ba
--- /dev/null
+++ b/scripts/ex2.watermassdistribution.jl
@@ -0,0 +1,39 @@
+#=
+ Example 2: Find the global distribution of a "water mass" defined by an oceanographically-relevant surface region.
+ Steps: (a) define the water mass 1). by a pre-defined surface
+            dyed with passive tracer concentration of 1,
+        (b) propagate the dye with the matrix A, with the result
+            being the fraction of water originating from the
+            surface region.
+ See Section 2b of Gebbie & Huybers 2010, esp. eqs. (15)-(17).
+=#
+using Revise
+using TMI, BenchmarkTools, PyPlot, PyCall
+
+pygui(true) #needed for Atom, not sure what it will do in other places
+
+TMIversion = "modern_180x90x33_GH10_GH12"
+TMIversion = "modern_90x45x33_GH10_GH12"
+A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)
+
+# choose water mass (i.e., surface patch) of interest
+# Enhancement: provide list of choices
+list = ("GLOBAL","ANT","SUBANT",
+            "NATL","NPAC","TROP","ARC",
+            "MED","ROSS","WED","LAB","GIN",
+            "ADEL","SUBANTATL","SUBANTPAC","SUBANTIND",
+            "TROPATL","TROPPAC","TROPIND")
+
+region = list[2]
+
+# do numerical analysis
+g = watermassdistribution(TMIversion,Alu,region,γ);
+
+# plot a section at 330 east longitude (i.e., 30 west)
+#lon_section = 329;
+lon_section = 330;
+gsection = section(g,lon_section,γ)
+lims = 0:5:100
+
+# make a plot of dye in the ocean
+dyeplot(γ.lat,-γ.depth[33:-1:1],100 * gsection[:,33:-1:1]', lims)
diff --git a/scripts/ex2.regeneration.jl b/scripts/ex3.regeneration.jl
similarity index 72%
rename from scripts/ex2.regeneration.jl
rename to scripts/ex3.regeneration.jl
index 5de7286f..e902cae3 100644
--- a/scripts/ex2.regeneration.jl
+++ b/scripts/ex3.regeneration.jl
@@ -9,12 +9,13 @@ pygui(true) #needed for Atom, not sure what it will do in other places
 
 # A, Alu, γ, inputfile = config(url,inputdir)
 # ΔPO₄ = readtracer(inputfile,"qpo4")
-TMIversion = "TMI_2010_2012_4x4x33"
-PO₄ᴿ, γ = regeneratedphosphate(TMIversion)
+TMIversion = "modern_90x45x33_GH10_GH12"
+A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)
+PO₄ᴿ = regeneratedphosphate(TMIversion,Alu,γ)
 
 
 # plot a section at 330 east longitude (i.e., 30 west)
-lon_section = 330;
+lon_section = 330; # only works if exact
 Psection = section(PO₄ᴿ,lon_section,γ)
 lims = 0:0.1:2.0
 
diff --git a/scripts/ex3.howoceansfilled.jl b/scripts/ex4.howoceansfilled.jl
similarity index 86%
rename from scripts/ex3.howoceansfilled.jl
rename to scripts/ex4.howoceansfilled.jl
index 4d5e2f45..1914774d 100644
--- a/scripts/ex3.howoceansfilled.jl
+++ b/scripts/ex4.howoceansfilled.jl
@@ -14,8 +14,10 @@
 using Revise
 using TMI, PyPlot, PyCall
 
-TMIversion = "TMI_2010_2012_4x4x33"
-volume = volumefilled(TMIversion)
+TMIversion = "modern_90x45x33_GH10_GH12"
+A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)
+
+volume = volumefilled(TMIversion,Alu,γ)
 
 # view the surface
 cntrs = 1:0.25:6
diff --git a/scripts/ex4.surfaceorigin.jl b/scripts/ex5.surfaceorigin.jl
similarity index 89%
rename from scripts/ex4.surfaceorigin.jl
rename to scripts/ex5.surfaceorigin.jl
index ea221392..c263c8ae 100644
--- a/scripts/ex4.surfaceorigin.jl
+++ b/scripts/ex5.surfaceorigin.jl
@@ -22,7 +22,10 @@ xlat = -6.38;  # deg N.
 xdepth = 3000;  # meters.
 loc = (xlon,xlat,xdepth)
 
-origin, γ = surfaceorigin(TMIversion,loc)
+TMIversion = "modern_90x45x33_GH10_GH12"
+A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)
+
+origin = surfaceorigin(loc, Alu, γ)
 
 # units: effective thickness in log10(meters)
 contourf(γ.lon,γ.lat,log10.(origin'))
diff --git a/scripts/ex5.filterdata.jl b/scripts/ex6.steadyclimatology.jl
similarity index 93%
rename from scripts/ex5.filterdata.jl
rename to scripts/ex6.steadyclimatology.jl
index fa8ccf42..fd6ddf95 100644
--- a/scripts/ex5.filterdata.jl
+++ b/scripts/ex6.steadyclimatology.jl
@@ -25,15 +25,15 @@ using Revise, TMI, GoogleDrive
 using PyPlot, PyCall, Test
 #, Distributions, LinearAlgebra,  Zygote, ForwardDiff, Optim
 
-TMIversion = "TMI_2010_2012_4x4x33"
-A, Alu, γ, inputfile = config(TMIversion) 
+TMIversion = "modern_90x45x33_GH10_GH12"
+A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)
 
 # first guess of change to surface boundary conditions
 # ocean values are 0
 u₀ = zeros(Float64,sum(γ.wet[:,:,1]))
 
 # take synthetic, noisy observations
-y, W⁻, ctrue = sample_observations(TMIversion,"θ")
+y, W⁻, ctrue = sample_observations(TMIversion,"θ",γ)
 
 # a first guess: observed surface boundary conditions are perfect.
 # set surface boundary condition to the observations.
@@ -44,7 +44,7 @@ d₀[:,:,1] = y[:,:,1]
 fg!(F,G,x) = costfunction_obs!(F,G,x,Alu,d₀,y,W⁻,γ)
 
 # filter the data with an Optim.jl method
-out = filterdata(u₀,Alu,y,d₀,W⁻,fg!,γ)
+out = steadyclimatology(u₀,Alu,y,d₀,W⁻,fg!,γ)
 
 # reconstruct by hand to double-check.
 ũ = out.minimizer
diff --git a/scripts/ex6.sparsedatamap.jl b/scripts/ex7.sparsedatamap.jl
similarity index 95%
rename from scripts/ex6.sparsedatamap.jl
rename to scripts/ex7.sparsedatamap.jl
index 46bc8688..c9e4d2e1 100644
--- a/scripts/ex6.sparsedatamap.jl
+++ b/scripts/ex7.sparsedatamap.jl
@@ -19,8 +19,8 @@
 =#
 using Revise, TMI, Interpolations
 
-TMIversion = "TMI_2010_2012_4x4x33"
-A, Alu, γ, TMIfile = config(TMIversion)
+TMIversion = "modern_90x45x33_GH10_GH12"
+A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)
 
 # first guess of change to surface boundary conditions
 # how many randomly sampled observations?
@@ -30,7 +30,7 @@ N = 20
 u₀ = zeros(Float64,sum(γ.wet[:,:,1]))
 
 # take synthetic, noisy observations
-y, W⁻, ctrue, locs, wis = sample_observations(TMIversion,"θ",N)
+y, W⁻, ctrue, locs, wis = sample_observations(TMIversion,"θ",γ,N)
 
 # make a silly first guess for surface
 d₀ = tracerinit(γ.wet)
diff --git a/scripts/mat2netcdf.jl b/scripts/mat2netcdf.jl
new file mode 100644
index 00000000..669e90fa
--- /dev/null
+++ b/scripts/mat2netcdf.jl
@@ -0,0 +1,26 @@
+# TMI solutions were originally saved in mat format. Here convert to NetCDF using 
+using Revise, MAT, LinearAlgebra, SparseArrays, TMI, DrWatson, NCDatasets, Test
+
+# choose one version
+TMIversion = "modern_90x45x33_GH10_GH12"
+TMIversion = "modern_180x90x33_GH10_GH12"
+TMIversion = "modern_90x45x33_unpub12"
+TMIversion = "modern_90x45x33_G14"
+
+# original NetCDF version
+#A, Alu, γ, TMIfile = config(TMIversion)
+
+@time A, Alu, γ, TMIfile, L, B = config_from_mat(TMIversion);
+
+config2nc(TMIversion,A,γ,L,B)
+
+filenetcdf = datadir("TMI_"*TMIversion*".nc")
+NCDataset(filenetcdf)
+
+# read from NetCDF?
+@time A2, Alu2, γ2, TMIfile2, L2, B2 = config_from_nc(TMIversion);
+
+# Do NetCDF and mat-files agree?
+@test A2 == A
+@test L2 == L
+@test B2 == B
diff --git a/src/TMI.jl b/src/TMI.jl
index df4e9afb..c0550ed0 100644
--- a/src/TMI.jl
+++ b/src/TMI.jl
@@ -4,27 +4,28 @@ using Revise
 using LinearAlgebra, SparseArrays, NetCDF, Downloads,
     GoogleDrive, Distances, DrWatson, GibbsSeaWater,  
     PyPlot, PyCall, Distributions, Optim,
-    Interpolations, LineSearches
+    Interpolations, LineSearches, MAT, NCDatasets
 
-export config, download,
-    vec2fld, fld2vec, depthindex, surfaceindex,
+export config, config_from_mat, config_from_nc,
+    vec2fld, fld2vec, surfaceindex,
+    lonindex, latindex, depthindex,
     surfacepatch, section,
     layerthickness, cellarea, cellvolume,
     planview, dyeplot, plotextent, tracerinit,
-    updateLinearindex,
-    watermassmatrixXYZ, watermassmatrixZYX,
-    linearindexXYZ, nearestneighbor,
+    watermassmatrix, watermassdistribution,
+    circulationmatrix, boundarymatrixXYZ,
+    linearindex, nearestneighbor, updatelinearindex,
     nearestneighbormask, horizontaldistance,
-    readtracer, cartesianindexZYX, Γ,
+    readtracer, cartesianindex, Γ,
     costfunction_obs, costfunction_obs!,
     costfunction, costfunction!,
     trackpathways, regeneratedphosphate, volumefilled,
-    surfaceorigin, sample_observations, filterdata,
+    surfaceorigin, sample_observations, steadyclimatology,
     steady_inversion,
     interpweights, interpindex,
     wetlocation, iswet,
     control2state, control2state!,
-    sparsedatamap
+    sparsedatamap, config2nc
 
 #Python packages - initialize them to null globally
 #const patch = PyNULL()
@@ -51,18 +52,18 @@ function __init__()
  end
 
 struct grid
-    lon::Vector{Real}
-    lat::Vector{Real}
-    depth::Vector{Real}
+    lon::Vector{Float64}
+    lat::Vector{Float64}
+    depth::Vector{Float64}
     I::Vector{CartesianIndex{3}} # index
-    R::Array{Integer,3}
-#        R::Array{Union{Integer,Nothing},3}
+    R::Array{Int,3}
 #    R::LinearIndices{3, Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}} 
     wet::BitArray{3}
 end
 
 """
-    function config(url,inputdir)
+    function config_from_nc(TMIversion)
+    Configure TMI environment from NetCDF input file.
 # Arguments
 - `TMIversion`: TMI version for water-mass/circulation model
 # Output
@@ -71,72 +72,143 @@ end
 - `γ`: TMI grid properties
 - `TMIfile`: TMI file name
 """
-function config(TMIversion)
+function config_from_nc(TMIversion)
 
     #- `url`: Google Drive URL for data
-    url = gdriveurl(TMIversion)
+    url = ncurl(TMIversion)
     
-    TMIfile = datadir(TMIversion*".nc")
+    TMIfile = datadir("TMI_"*TMIversion*".nc")
+    println(url)
+    println(TMIfile)
     !isdir(datadir()) ? mkpath(datadir()) : nothing
-    !isfile(TMIfile) ? download(url,datadir()) : nothing
+    !isfile(TMIfile) ? google_download(url,datadir()) : nothing
 
-    ncdata = NetCDF.open(TMIfile)
-    #println(ncdata)
-    
-    # move this to runtests.jl to see if it is read correctly
-    # Azyx = watermassmatrixZYX(TMIfile)
+    #ncdata = NetCDF.open(TMIfile) # necessary?
 
-    # make a sample field from zyx cartesian indices
-    Izyx = cartesianindexZYX(TMIfile)
+    # read Cartesian Index from file.
+    I = cartesianindex(TMIfile)
 
     # make a mask
+    wet = falses(maximum(I)[1],maximum(I)[2],maximum(I)[3])
+    #wet = BitArray{3}(undef,maximum(I)[1],maximum(I)[2],maximum(I)[3])
+    #fill!(wet,0)
+    wet[I] .= 1
+
+    R = linearindex(wet)
+
+    println("A")
+    @time A = watermassmatrix(TMIfile)
+
+    # LU factorization for efficient matrix inversions
+    println("Alu")
+    @time Alu = lu(A)
+    
+    # get properties of grid
+    lat,lon,depth = gridprops(TMIfile)
+    γ = grid(lon,lat,depth,I,R,wet)
+
+    # would be good to make this optional
+    println("L=")
+    @time L = circulationmatrix(TMIfile,A,γ)
+
+    println("B=")
+    @time B = boundarymatrix(TMIfile,γ)
+    
+    return  A, Alu, γ, TMIfile, L, B
+
+end
+
+"""
+Configure TMI environment from original MATLAB output
+"""
+function config_from_mat(TMIversion)
+
+    #- `url`: Google Drive URL for data
+    url = maturl(TMIversion)
+    TMIfile = datadir("TMI_"*TMIversion*".mat")
+    TMIfilegz = TMIfile*".gz"
+    println(TMIfile)
+    !isdir(datadir()) ? mkpath(datadir()) : nothing
+    !isfile(TMIfilegz) & !isfile(TMIfile) ? google_download(url,datadir()) : nothing
+
+    # cloak mat file in gz to get Google Drive spam filter to shut down
+    isfile(TMIfilegz) & !isfile(TMIfile) ? run(`gunzip $TMIfilegz`) : nothing
+    
+    # move this to runtests.jl to see if it is read correctly?
+    # Azyx = watermassmatrix(TMIfile) 
+
+    # # make a sample field from zyx cartesian indices
+    Izyx = cartesianindex(TMIfile)
+
+    # # make a mask
     wet = BitArray{3}(undef,maximum(Izyx)[1],maximum(Izyx)[2],maximum(Izyx)[3])
     fill!(wet,0)
     wet[Izyx] .= 1
 
-    # if a tracer is available, should be consistent with this definition
-    #wet = .!isnan.(c)
-    
-    # need to write this function
-    I = cartesianindexXYZ(wet)
+    # # consistent with tracer definition?
+    # #wet = .!isnan.(c)
 
-    R = linearindexXYZ(wet)
+    I = cartesianindex(wet)
 
-    A = watermassmatrixXYZ(TMIfile,R)
-    #R = R[ wet ] # eliminate land points
+    R = linearindex(wet)
 
-    # LU factorization for efficient matrix inversions
+    Azyx = watermassmatrix(TMIfile)
+    A = Azyx2xyz(TMIfile,Azyx,R)
+    # #R = R[ wet ] # eliminate land points
+
+    # # LU factorization for efficient matrix inversions
     Alu = lu(A)
-    
+
     # get properties of grid
     lat,lon,depth = gridprops(TMIfile)
 
     γ = grid(lon,lat,depth,I,R,wet)
 
-    return  A, Alu, γ, TMIfile
-
+    # need to make this optional
+    L = circulationmatrix(TMIfile,γ)
+    
+    B = boundarymatrix(TMIfile,γ)
+    
+    # consider re-ordering this.
+    # some output should be optional
+    # return Izyx or I or neither?
+    #return  A, Alu, γ, TMIfile, I, L, B
+    return  A, Alu, γ, TMIfile, L, B
 end
 
 """
-    function cartesianindexZYX(file)
+    function cartesianindex(file)
     Read and assemble the grid coordinates
     according to the legacy MATLAB code (z,y,x order).
 # Arguments
 - `file`: TMI NetCDF file name
 # Output
-- `grid`: TMI grid coordinates
+- `I`: TMI Cartesian index for wet points
 """
-function cartesianindexZYX(file)
+function cartesianindex(file::String)
     # make the Cartesian tracer grid
-    it = convert(Array{Int,1},ncread(file,"xgrid"))
-    jt = convert(Array{Int,1},ncread(file,"ygrid"))
-    kt = convert(Array{Int,1},ncread(file,"zgrid"))
-    I = CartesianIndex.(it,jt,kt)
+    if file[end-1:end] == "nc"
+
+        it = convert(Vector{Int},ncread(file,"i"))
+        jt = convert(Vector{Int},ncread(file,"j"))
+        kt = convert(Vector{Int},ncread(file,"k"))
+
+        I = CartesianIndex.(it,jt,kt)
+
+    elseif file[end-2:end] == "mat"
+
+        matobj = matopen(file)
+        haskey(matobj,"it") ? it=convert(Vector{Integer},vec(read(matobj,"it"))) : it = convert(Vector{Integer},vec(read(matobj,"i")))
+        haskey(matobj,"jt") ? jt=convert(Vector{Integer},vec(read(matobj,"jt"))) : jt=convert(Vector{Integer},vec(read(matobj,"j")))
+        haskey(matobj,"kt") ? kt=convert(Vector{Integer},vec(read(matobj,"kt"))) : kt=convert(Vector{Integer},vec(read(matobj,"k")))
+        close(matobj)
+        I = CartesianIndex.(it,jt,kt) # should this be reversed?
+    end
     return I
 end
 
 """
-    function cartesianindexXYZ(wet)
+    function cartesianindex(wet)
     Read and assemble the grid coordinates
     according to a 3D tracer in x,y,z order
 # Arguments
@@ -144,17 +216,17 @@ end
 # Output
 - `I`: 3D Cartesian indices
 """
-cartesianindexXYZ(wet) = findall(wet)
+cartesianindex(wet::BitArray{3}) = findall(wet)
 
 """
-    function linearindexXYZ(file)
+    function linearindex(wet)
     Read and assemble the grid coordinates.
 # Arguments
 - `wet`: 3D mask for wet points
 # Output
 - `R`: array of linear indices, but not a LinearIndices type
 """
-function linearindexXYZ(wet)
+function linearindex(wet)
     R = Array{Int64,3}(undef,size(wet))
     fill!(R,0)
     # R = Array{Union{Int64,Nothing},3}(nothing,size(wet))
@@ -174,183 +246,242 @@ end
 - `grid`: TMI grid coordinates
 """
 function gridprops(file)
-    lat = ncread(file,"lat")
-    lon = ncread(file,"lon")
-    depth = ncread(file,"depth")
+    if file[end-1:end] == "nc"
+        
+        lat = convert(Vector{Float64},ncread(file,"lat"))
+        lon = convert(Vector{Float64},ncread(file,"lon"))
+        depth = convert(Vector{Float64},ncread(file,"depth"))
+
+    elseif file[end-2:end] == "mat"
+        
+        matobj = matopen(file)
+        lon=convert(Vector{Float64},vec(read(matobj,"LON")))
+        lat=convert(Vector{Float64},vec(read(matobj,"LAT")))
+        depth=convert(Vector{Float64},vec(read(matobj,"DEPTH")))
+        close(matobj)
+
+    end
+    
     return lat,lon,depth
 end
 
 """
-    function watermassmatrixZYX(file)
+    function watermassmatrix(file)
     Read and assemble the water-mass matrix.
-    Legacy version from MATLAB.
 # Arguments
-- `file`: TMI NetCDF file name
+- `file`: TMI NetCDF or MATLAB file name
 # Output
 - `A`: water-mass matrix
 """
-function watermassmatrixZYX(file)
-    i = ncread(file,"i")
-    j = ncread(file,"j")
-    m = ncread(file,"m")
-    A = sparse(i,j,m)
+function watermassmatrix(file)
+
+    # consider adding a catch if A doesn't exist in file.
+
+    if file[end-1:end] == "nc"
+        # Int or Integer?
+        i = convert(Vector{Int},ncread(file,"Arow"))
+        j = convert(Vector{Int},ncread(file,"Acol"))
+        m = ncread(file,"m")
+        A = sparse(i,j,m)
+    elseif file[end-2:end] == "mat"
+        matobj = matopen(file)
+        if haskey(matobj,"A")
+            A=read(matobj,"A")
+            close(matobj)
+        else
+            close(matobj)
+            return nothing
+        end
+
+        # But MATLAB had zyx format and we need xyz format.
+        # linearindices R not available so will do conversion in higher scope
+
+    end
     return A
 end
 
 """
-        function watermassmatrixXYZ(file,R)
-    Read and assemble the water-mass matrix from MATLAB.
-    Transfer to updated x,y,z version
+        function Azyx2xyz(TMIfile,Azyx,γ)
+   
+    Transfer zyx format water-mass matrix A to xyz format
 # Arguments
-- `file`: TMI NetCDF file name
+- `Azyx`: water-mass matrix in zyx format
 - `γ`: TMI grid
 # Output
-- `A`: water-mass matrix
+- `Axyz`: water-mass matrix in xyz format
 """
-function watermassmatrixXYZ(file,R)
-
-    # MATLAB accounting z,y,x
-    izyx = convert(Vector{Int},ncread(file,"i"))
-    jzyx = convert(Vector{Int},ncread(file,"j"))
+function Azyx2xyz(file,Azyx,R)
 
+    izyx, jzyx, mzyx = findnz(Azyx)
+    Izyx = cartesianindex(file)
+        
     # Julia accounting x,y,z
-    Izyx = cartesianindexZYX(file)
     ixyz = updatelinearindex(izyx,Izyx,R)
     jxyz = updatelinearindex(jzyx,Izyx,R)
     
-    m = ncread(file,"m")
-
     # use grid indices to switch i,j values
-    A = sparse(ixyz,jxyz,m)
-    #A = A[wet,:];
-    return A
+    Axyz = sparse(ixyz,jxyz,mzyx)
+    return Axyz
 end
 
-
 """
-    function updatelinearindex(izyx,Izyx,R)
-    Linear index translated from z,y,x to x,y,z accounting
+    function circulationmatrix(file,γ)
+    Read and assemble the circulation matrix from MATLAB.
+    Transfer to updated x,y,z version
 # Arguments
-- `izyx`: index of interest in z,y,x accounting
-- `Izyx`: wet Cartesian Index for z,y,x
-- `R`: Linear indices for x,y,z 
+- `file`: TMI MATLAB file name
+- `γ`: TMI grid
 # Output
-- `ixyz`: index of interest in x,y,z accounting
+- `L`: circulation matrix in xyz format
 """
-function updatelinearindex(izyx,Izyx,R)
-    # get Izyx Cartesian index stored from legacy MATLAB code
-    ixyz = R[Izyx[izyx]]
-    return ixyz
+function circulationmatrix(file,γ)
+
+    if file[end-2:end] == "mat" 
+
+        matobj = matopen(file)
+        if haskey(matobj,"L")
+            # Matlab output in zyx format
+            Lzyx=read(matobj,"L")
+            close(matobj)
+
+            Izyx = cartesianindex(file)
+            izyx, jzyx, Fzyx = findnz(Lzyx)
+            # Julia accounting x,y,z
+            ixyz = updatelinearindex(izyx,Izyx,γ.R)
+            jxyz = updatelinearindex(jzyx,Izyx,γ.R)
+            L = sparse(ixyz,jxyz,Fzyx)
+
+        else
+            close(matobj)
+            return nothing
+        end
+
+    elseif file[end-1:end] == "nc"
+
+        # based on function arguments, read from inefficient storage of L matrix.
+        if haskey(NCDataset(file),"F")
+
+            i = convert(Vector{Int},ncread(file,"Lrow"))
+            j = convert(Vector{Int},ncread(file,"Lcol"))
+            F = ncread(file,"F")
+            L = sparse(i,j,F)
+        else
+            return nothing
+        end
+    end
+    
+    return L
 end
 
 """
-    function watermassmatrix(file)
-    Read and assemble the water-mass matrix.
+    function circulationmatrix(file,A,γ)
+    Read and assemble the circulation matrix from the efficient storage of A and F₀ variables. 
 # Arguments
-- `file`: TMI NetCDF file name
+- `file`: TMI MATLAB file name
+- `A`: TMI water-mass matrix
+- `γ`: TMI grid
 # Output
-- `A`: water-mass matrix
+- `L`: circulation matrix in xyz format
 """
-function watermassmatrix(file)
+function circulationmatrix(file,A,γ)
 
-    m = massFractions(file)
+    file[end-1:end] !== "nc" && error("not a NetCDF file")
 
-    A = watermassmatrix(m)
-    
-    # assemble m into A
-    # Atest = Array{SparseArrays.SparseMatrixCSC{Float64, Int64},3}
+    # based on function arguments, read F₀ to efficiently reproduce L matrix.
 
+    if haskey(NCDataset(file),"F₀")
+        F₀ = ncread(file,"F₀")
+        F₀vec = F₀[γ.wet]
     
-    # it = convert(Vector{Int},ncread(file,"xgrid"))
-    # jt = convert(Vector{Int},ncread(file,"ygrid"))
-    # kt = convert(Vector{Int},ncread(file,"zgrid"))
-    # i = convert(Vector{Int},ncread(file,"i"))
-    # j = convert(Vector{Int},ncread(file,"j"))
-    # m = ncread(file,"m")
+        # For each row of A, multiply by F₀
+        i, j, F = findnz(A)
 
-    # for gg = 1:length(i)
-
-    #     list = findall(x -> x == gg, i)
+        # careful, this loop can be really slow
+        for nn in eachindex(i)
+            F[nn] *= F₀vec[i[nn]]
+        end
 
-    #     # get linear coordinates.
-    #     [linear[q] = R[it[list[q]],jt[list[q]],kt[list[q]]] for q = 1:length(list)] 
-        
-    #     # can't store sparse 3D matrix. just 2D.
-    #     # need to go from i,j,k to linear index
-    #     Arowsparse = sparse(  m[gg])
-    #     A[it[i[gg]],jt[i[gg]],kt[i[gg]]] = sparse([it[j[gg]],jt[j[gg]],kt[j[gg]]] = m[gg]
-    # end
+        L = sparse(i,j,F)
+        return L
+    else
+        return nothing
+    end
     
-    #A = sparse(i,j,m)
-    return A
 end
 
 """
-    function watermassmatrix(file)
-    Assemble the water-mass matrix given mass fractions `m`
+        function boundarymatrix(file,γ)
+    Read and assemble the boundary matrix from MATLAB.
+    Transfer to updated x,y,z version
 # Arguments
-- `m`: mass fractions
+- `file`: TMI MATLAB file name
+- `γ`: TMI grid
 # Output
-- `A`: water-mass matrix
-"""
-# function watermassmatrix(m)
-#     # Goal:assemble m into A
-
-#     # preallocate A
-#     A = Array{SparseArrays.SparseMatrixCSC{Float64, Int64},3}
-                                                  
-#     # loop over each equation
-    
-#     # get linear index for equation (destination)
+- `B`: boundary condition matrix
+"""
+function boundarymatrix(file,γ)
+
+    if file[end-2:end] == "mat"
+
+        matobj = matopen(file)
+        if haskey(matobj,"B")
+            Bzyx=read(matobj,"B")
+            close(matobj)
+
+            # matlab in zyx format.
+            # consider using Azyx2xyz here.
+            Izyx = cartesianindex(file)
+            izyx, jzyx, Fzyx = findnz(Bzyx)
+            # for B, rows are 3D grid space, columns are for the surface index. 
+            # Julia accounting x,y,z
+            Isfc = surfaceindex(Izyx)
+            ixyz = updatelinearindex(izyx,Izyx,γ.R)
+            jxyz = updatelinearindex(Isfc[jzyx],Izyx,γ.R)
+
+            # assume surface at k = 1 (revisit for LGM problem)
+            # give the full dimension of sparse matrix
+            B = sparse(ixyz,jxyz,Fzyx,sum(γ.wet),sum(γ.wet[:,:,1]))
+        else
+            close(matobj)
+            return nothing
+        end
 
-#     # get linear indices for sources
+    elseif file[end-1:end] == "nc"
 
-#     # make a list of row (destination), column (destination), m value
+        # based on function arguments, read from inefficient storage of L matrix.
+        if haskey(NCDataset(file),"b")
 
-#     # complete loop
+            i = convert(Vector{Int},ncread(file,"Brow"))
+            j = convert(Vector{Int},ncread(file,"Bcol"))
+            b = ncread(file,"b")
+            B = sparse(i,j,b,sum(γ.wet),sum(γ.wet[:,:,1]))
 
-#     # make sparse matrix
-#     #A = sparse(i,j,m)
-#     return A
-# end
+        else
+            return nothing
+        end
+    end
+    return B
+end
 
 """
-    function massFractions(file)
-    Read and assemble the water-mass fractions, `m`
+    function updatelinearindex(izyx,Izyx,R)
+    Linear index translated from z,y,x to x,y,z accounting
 # Arguments
-- `file`: TMI NetCDF file name
+- `izyx`: index of interest in z,y,x accounting
+- `Izyx`: wet Cartesian Index for z,y,x
+- `R`: Linear indices for x,y,z 
 # Output
-- `m`: mass fractions
+- `ixyz`: index of interest in x,y,z accounting
 """
-# function massFractions(file)
-#     Atest = Array{SparseArrays.SparseMatrixCSC{Float64, Int64},3}
-    
-#     it = convert(Vector{Int},ncread(file,"xgrid"))
-#     jt = convert(Vector{Int},ncread(file,"ygrid"))
-#     kt = convert(Vector{Int},ncread(file,"zgrid"))
-#     i = convert(Vector{Int},ncread(file,"i"))
-#     j = convert(Vector{Int},ncread(file,"j"))
-#     mMat = ncread(file,"m") # MATLAB generated list
-
-#     for gg = 1:length(i)
-
-#         list = findall(x -> x == gg, i)
-
-#         # get linear coordinates.
-#         [linear[q] = R[it[list[q]],jt[list[q]],kt[list[q]]] for q = 1:length(list)] 
-        
-#         # can't store sparse 3D matrix. just 2D.
-#         # need to go from i,j,k to linear index
-#         Arowsparse = sparse(  m[gg])
-#         A[it[i[gg]],jt[i[gg]],kt[i[gg]]] = sparse([it[j[gg]],jt[j[gg]],kt[j[gg]]] = m[gg]
-#     end
-    
-#     return m
-# end
+function updatelinearindex(izyx,Izyx,R)
+    # get Izyx Cartesian index stored from legacy MATLAB code
+    ixyz = R[Izyx[izyx]]
+    return ixyz
+end
                                                   
 """
     function readtracer(file,tracername)
-    Read and assemble the water-mass matrix.
+    Read a tracer field from NetCDF.
 # Arguments
 - `file`: TMI NetCDF file name
 - `tracername`: name of tracer
@@ -362,6 +493,9 @@ function readtracer(file,tracername)
     return c
 end
 
+"""
+    Horizontal area of grid cell
+"""
 function cellarea(γ)
     dx = zonalgriddist(γ)
     dy = haversine((γ.lon[1],γ.lat[1])
@@ -380,6 +514,9 @@ function cellarea(γ)
     return area
 end
 
+"""
+    Volume of each grid cell.
+"""
 function cellvolume(γ)
     dz = layerthickness(γ)
     area = cellarea(γ)
@@ -410,28 +547,6 @@ function zonalgriddist(γ::grid)
     return dx
 end
 
-"""
-    function download(url,inputdir)
-    Read and assemble all TMI inputs.
-# Arguments
-- `url`: Google Drive location of TMI input
-- `inputdir`: input directory location to store file
-# Output
-- none
-"""
-function download(url,inputdir)
-    # later include options and move these settings to arguments.
-
-    # make sure input dir exists
-    !isdir(inputdir) ? mkdir(inputdir) : nothing
-
-    # two ways to download
-    # 1. Use `run` for a shell command (less portable). Also difficult for Google Drive. See downloadTMIfromGoogleDrive.sh.
-
-    # 2. Use GoogleDrive.jl package
-    google_download(url,inputdir)
-end
-
 """
     function vec2fld
     Transfer a vector to a 3D field with accounting for ocean bathymetry
@@ -441,12 +556,17 @@ end
 # Output
 - `field`: field in 3d form including land points (NaN)
 """
-function vec2fld(vector::Vector{Float64},I::Vector{CartesianIndex{3}})
+function vec2fld(vector::Vector{T},I::Vector{CartesianIndex{3}}) where T<:Real
+
+    # choose NaN for now, zero better? nothing better?
+    fillvalue = zero(T)/zero(T) # NaN32 or NaN64
 
     nx = maximum(I)[1]
     ny = maximum(I)[2]
     nz = maximum(I)[3]
-    field = NaN .* zeros(nx,ny,nz)
+
+    # faster instead to allocate as undef and then fill! ?
+    field = fillvalue .* zeros(nx,ny,nz)
 
     # a comprehension
     [field[I[n]]=vector[n] for n ∈ eachindex(I)]
@@ -455,15 +575,17 @@ end
 
 """
     function fld2vec
-    Transfer 3D field with accounting for ocean bathymetry to a vector without land points
+    Transfer 3D field with accounting for ocean bathymetry to a vector without land points.
+    This is done more easily with a BitArray mask, i.e., vector = field[mask].
+    This function may be removed in the future.
 # Arguments
 - `field`: field in 3d form including land points (NaN)
 - `I`: cartesian indices of ocean points
 # Output
 - `vector`: field in vector form (no land points)
 """
-function fld2vec(field::Array{Float64,3},I::Vector{CartesianIndex{3}})
-    vector = Vector{Real}(undef,length(I))
+function fld2vec(field::Array{T,3},I::Vector{CartesianIndex{3}}) where T<:Real
+    vector = Vector{T}(undef,length(I))
     #- a comprehension
      [vector[n] = field[I[n]] for n ∈ eachindex(I)];
      return vector
@@ -480,7 +602,7 @@ function fld2vec(field::Array{Float64,3},I::Vector{CartesianIndex{3}})
 # Output
 - `d`: vector that describes surface patch
 """
-function surfacepatch(lonbox::Vector{T},latbox::Vector{T},γ::grid)::Array{Float64} where T<:Real
+function surfacepatch(lonbox,latbox,γ::grid)
 
     # ternary operator to handle longitudinal wraparound
     lonbox[1] ≤ 0 ? lonbox[1] += 360 : nothing
@@ -489,10 +611,10 @@ function surfacepatch(lonbox::Vector{T},latbox::Vector{T},γ::grid)::Array{Float
     # define the surface boundary condition
 
     # preallocate
-    d = tracerinit(γ.wet)
+    d = tracerinit(γ.wet,Float64)
 
+    # can you add a logical to a Float64? yes, it's 1.0
     [d[i,j,1] =  latbox[1] ≤ γ.lat[j] ≤ latbox[2] && lonbox[1] ≤ γ.lon[i] ≤ lonbox[2] for i in eachindex(γ.lon) for j in eachindex(γ.lat)] 
-    d[.!γ.wet] .= NaN # double check that NaNs stay NaNs
 
     # old method for vectors
         #nfield = length(γ.I) # number of ocean points
@@ -576,9 +698,13 @@ end
 """
 function horizontaldistance(loc,γ::grid)
 
+    # hordist will have same type as lon,lat,depth
+    T = eltype(γ.lon)
+    
     # pre-allocate horizontal distance
-    hordist = Matrix{Float64}(undef,length(γ.lon),length(γ.lat))
-    fill!(hordist,NaN)
+    hordist = Matrix{T}(undef,length(γ.lon),length(γ.lat))
+    # will give NaN with appropriate precision
+    fill!(hordist,zero(T)/zero(T))
     
     # calculate haversine horizontal distance on sphere
     [hordist[γ.I[ii]] = haversine((loc[1],loc[2]),                  (γ.lon[γ.I[ii][1]],γ.lat[γ.I[ii][2]]))
@@ -683,12 +809,35 @@ end
     function depthindex(I) 
     Get the k-index (depth level) from the Cartesian index
 """
-function depthindex(I) 
-    k = Vector{Int64}(undef,length(I))
-    [k[n]=I[n][3] for n ∈ 1:length(I)]
+function depthindex(I)
+    T = eltype(I[1])
+    k = Vector{T}(undef,length(I))
+    [k[n]=I[n][3] for n ∈ eachindex(I)]
     return k
 end
 
+"""
+    function lonindex(I) 
+    Get the i-index (lon index) from the Cartesian index
+"""
+function lonindex(I)
+    T = eltype(I[1])
+    i = Vector{T}(undef,length(I))
+    [i[n]=I[n][1] for n ∈ eachindex(I)]
+    return i
+end
+
+"""
+    function latindex(I) 
+    Get the j-index (latitude index) from the Cartesian index
+"""
+function latindex(I)
+    T = eltype(I[1])
+    j = Vector{T}(undef,length(I))
+    [j[n]=I[n][2] for n ∈ eachindex(I)]
+    return j
+end
+
 """
     function surfaceindex(I) 
     Get the vector-index where depth level == 1 and it is ocean.
@@ -704,6 +853,7 @@ end
     perhaps better to have a tracer struct and constructor
 # Arguments
 - `wet`::BitArray mask of ocean points
+- `ltype`:: optional type argument, default=Float64
 # Output
 - `d`:: 3d tracer field with NaN on dry points
 """
@@ -718,6 +868,29 @@ function tracerinit(wet,ltype=Float64)
     return d
 end
 
+""" 
+    function tracerinit(wet,vec,I)
+          initialize tracer field on TMI grid
+        perhaps better to have a tracer struct and constructor
+# Arguments
+- `wet`:: BitArray mask of ocean points
+- `vec`:: vector of values at wet points
+- `I`:: Cartesian Index for vector
+# Output
+- `field`:: 3d tracer field with NaN on dry points
+"""
+function tracerinit(vec,I,wet)
+
+    # preallocate
+    T = eltype(vec)
+    field = Array{T}(undef,size(wet))
+    fill!(field,zero(T)/zero(T))    
+
+    #- a comprehension
+    [field[I[n]]=vec[n] for n ∈ eachindex(I)]
+    return field
+end
+
 """ 
     function control2state(tracer2D,γ)
     turn 2D surface field into 3D field with zeroes below surface    
@@ -793,6 +966,8 @@ function state2obs(cvec,wis,γ)
     N = length(wis)
     sumwis = Vector{Float64}(undef,N)
     list = vcat(1:length(γ.lon),1)
+
+    # perhaps the most clever line in TMI.jl?
     wetwrap = view(γ.wet,list,:,:)
 
     [sumwis[i] = wetwrap[wis[i]...] for i in eachindex(wis)]
@@ -802,13 +977,11 @@ function state2obs(cvec,wis,γ)
     c̃ = tracerinit(γ.wet)
     c̃[γ.wet] = cvec
     replace!(c̃,NaN=>0.0)
-    [ỹ[i] = c̃[wis[i]...]/sumwis[i] for i in 1:N]
+    cwrap = view(c̃,list,:,:)
+
+    [ỹ[i] = cwrap[wis[i]...]/sumwis[i] for i in eachindex(wis)]
     return ỹ
 end
-
-# can I use multiple dispatch like this?
-# don't think so.
-#Γ(x) = Γ(x,wet)
     
 """ 
     function trackpathways(TMIversion,latbox,lonbox)
@@ -819,13 +992,11 @@ end
 - `TMIversion`: version of TMI water-mass/circulation model
 - `latbox`: min and max latitude of box
 - `lonbox`: min and max longitude of box
+- `γ`: TMI grid
 # Output
 - `c`: fraction of water from surface source
-- `γ`: TMI grid
 """
-function trackpathways(TMIversion,latbox,lonbox)
-
-    A, Alu, γ = config(TMIversion)
+function trackpathways(Alu,latbox,lonbox,γ)
 
     d = surfacepatch(lonbox,latbox,γ)
 
@@ -833,47 +1004,107 @@ function trackpathways(TMIversion,latbox,lonbox)
     c = tracerinit(γ.wet); # pre-allocate c
 
     # make methods that make the "wet" index unnecessary
-    c[γ.wet] = Alu\d[γ.wet] # presumably equivalent but faster than `c = A\d`
+    c[γ.wet] = Alu\d[γ.wet] # equivalent but faster than `c = A\d`
 
-    return c, γ
+    return c
 end
 
 """ 
-    function gdriveurl(TMIversion)
-    placeholder function to give location (URL) of Google Drive input
+    function watermassdistribution(TMIversion,latbox,lonbox)
+    Track the pathways of a user-defined water mass.
+     Steps: (a) define the water mass by an oceanographically-relevant surface patch dyed with passive tracer concentration of one
+         (b) propagate the dye with the matrix A, with the result being the fraction of water originating from the surface region.
+     See Section 2b of Gebbie & Huybers 2010, esp. eqs. (15)-(17).
+# Arguments
+- `TMIversion`: version of TMI water-mass/circulation model
+- `Alu`: LU decomposition of water-mass matrix A
+- `region`: name of pre-defined surface region
+- `γ`: TMI grid
+# Output
+- `g`: water-mass fraction
+"""
+function watermassdistribution(TMIversion,Alu,region,γ)
+
+    d = surfaceregion(TMIversion,region,γ)
+
+    # do matrix inversion to get quantity of dyed water throughout ocean:
+    g = tracerinit(γ.wet); # pre-allocate c
+
+    # make methods that make the "wet" index unnecessary
+    g[γ.wet] = Alu\d[γ.wet] # equivalent but faster than `c = A\d`
+
+    return g
+end
+
+""" 
+    function ncurl(TMIversion)
+    placeholder function to give location (URL) of NetCDF Google Drive input
     in the future, consider a struct or Dict that describes all TMI versions.
 # Arguments
 - `TMIversion`: version of TMI water-mass/circulation model
 # Output
 - `url`: location (URL) for download
 """
-function gdriveurl(TMIname)
-    if TMIname == "TMI_2010_2012_4x4x33"
-        #        url = "https://docs.google.com/uc?export=download&id=1Zycnx6_nifRrJo8XWMdlCFv4ODBpi-i7"
+function ncurl(TMIname)
+    if TMIname == "modern_90x45x33_GH10_GH12"
         url = "https://docs.google.com/uc?export=download&id=1Fn_cY-90_RDbBGh6kV0kpXmsvwdjp1Cd"
+    elseif TMIname == "modern_180x90x33_GH10_GH12"
+        url = "https://docs.google.com/uc?export=download&id=1-YEkB_YeQGqPRH6kauhBb2bi_BjVGt9b"
+    elseif TMIname == "modern_90x45x33_unpub12"
+        url = "https://docs.google.com/uc?export=download&id=1Kw_Mr7fiKqan0nx0dKvGHnSInP0hQ7AV"
+    elseif TMIname == "modern_90x45x33_G14"
+        url = "https://docs.google.com/uc?export=download&id=1aeE7EXA-vy3Cm_drt4qCFw4AlpYrdudk"
     else
         url = nothing
     end
 end
 
 """ 
-    function regeneratedphosphate(TMIversion)
+    function maturl(TMIversion)
+    Find *mat file here.
+    placeholder function to give location (URL) of Google Drive input
+    in the future, consider a struct or Dict that describes all TMI versions.
+# Arguments
+- `TMIversion`: version of TMI water-mass/circulation model
+# Output
+- `url`: location (URL) for download
+"""
+function maturl(TMIname)
+    if TMIname == "modern_90x45x33_GH10_GH12"
+        url = "https://docs.google.com/uc?export=download&id=1qPRq7sonwdjkPhpMcAnvuv67OMZSBqgh"
+    elseif TMIname == "modern_180x90x33_GH10_GH12"
+        url = "https://docs.google.com/uc?export=download&id=11zD1nOfT6V7G0qIHdjK2pDGHFk-ExXwU"
+    elseif TMIname == "modern_90x45x33_unpub12"
+        url = "https://docs.google.com/uc?export=download&id=1sqkjFCPxZT_2Bm9rsp0acyxxkBri9YAT"
+    elseif TMIname == "modern_90x45x33_G14"
+        url = "https://docs.google.com/uc?export=download&id=1dCrDe5VXrsXiOf04mbuHID7xc5Ymm8-z"
+    else
+        url = nothing
+    end
+    return url
+end
+
+""" 
+    function regeneratedphosphate(TMIversion,Alu,γ)
     Regenerated (i.e., accumulated, remineralized) phosphate
 # Arguments
 - `TMIversion`: version of TMI water-mass/circulation model
+- `Alu`: LU decomposition of water-mass matrix A
+- `γ`: TMI grid
 # Output
 - `PO₄ᴿ`: regenerated phosphate
-- `γ`: TMI grid
 """
-function regeneratedphosphate(TMIversion)
+function regeneratedphosphate(TMIversion,Alu,γ)
 
-    A, Alu, γ, inputfile = config(TMIversion)
-    ΔPO₄ = readtracer(inputfile,"qpo4")
+    inputfile = datadir("TMI_"*TMIversion*".nc")
+        
+    #A, Alu, γ, inputfile = config(TMIversion)
+    qPO₄ = readtracer(inputfile,"qPO₄")
 
     # PO₄ᴿ = cumulative regenerated phosphate
     PO₄ᴿ = tracerinit(γ.wet); # pre-allocate 
-    PO₄ᴿ[γ.wet] = -(Alu\ΔPO₄[γ.wet])
-    return PO₄ᴿ, γ
+    PO₄ᴿ[γ.wet] = -(Alu\qPO₄[γ.wet])
+    return PO₄ᴿ
 end
 
 """ 
@@ -888,12 +1119,14 @@ end
      See Section 3 and Supplementary Section 4, Gebbie & Huybers 2011. 
 # Arguments
 - `TMIversion`: version of TMI water-mass/circulation model
+- `Alu`: LU decomposition of water-mass matrix A
+- `γ`: TMI.grid
 # Output
 - `volume`: global ocean volume filled by a surface region
 """
-function volumefilled(TMIversion)
+function volumefilled(TMIversion,Alu,γ)
 
-    A, Alu, γ = config(TMIversion)
+    #A, Alu, γ = config(TMIversion)
     
     v = cellvolume(γ)
     area = cellarea(γ)
@@ -904,8 +1137,9 @@ function volumefilled(TMIversion)
 
     # scale the sensitivity value by surface area so that converging meridians are taken into account.
     I = γ.I
-    volume = Matrix{Float64}(undef,length(γ.lon),length(γ.lat))
-    fill!(volume,0.0)
+    volume = zeros(Float64,length(γ.lon),length(γ.lat))
+    #volume = Matrix{Float64}(undef,length(γ.lon),length(γ.lat))
+    #fill!(volume,0.0)
 
     # this step could use a function with γ.I argument
     [volume[I[ii][1],I[ii][2]] = dVdd[I[ii]] / area[I[ii][1],I[ii][2]] for ii ∈ eachindex(I) if I[ii][3] == 1]
@@ -924,15 +1158,15 @@ end
     The sensitivity is exactly the mass fraction originating from each source.      
     This problem is mathematically similar to determining how the ocean is filled.
 # Arguments
-- `TMIversion`: version of TMI water-mass/circulation model
 - `loc`: location (lon,lat,depth) of location of interest
+- `Alu`: LU decomposition of water-mass matrix A
+- `γ`: TMI grid
 # Output
 - `origin`: surface map of fraction of source water for a given location
-- `γ`: TMI grid
 """
-function surfaceorigin(TMIversion,loc)
+function surfaceorigin(loc,Alu,γ)
 
-    A, Alu, γ = config(TMIversion)
+    #A, Alu, γ = config(TMIversion)
 
     ctmp = tracerinit(γ.wet)
     δ = interpweights(loc,γ)
@@ -947,7 +1181,7 @@ function surfaceorigin(TMIversion,loc)
 
     # origin is defined at sea surface
     origin = view(dvlocdd,:,:,1)
-    return origin, γ
+    return origin
 end
 
 """
@@ -1028,7 +1262,7 @@ function interpweights(loc,γ)
 end
 
 """
-function filterdata(u₀,Alu,y,d₀,W⁻,γ)
+function steadyclimatology(u₀,Alu,y,d₀,W⁻,γ)
      Find the distribution of a tracer given:
      (a) the pathways described by A or its LU decomposition Alu,
      (b) first-guess boundary conditions and interior sources given by d₀,
@@ -1048,7 +1282,7 @@ function filterdata(u₀,Alu,y,d₀,W⁻,γ)
 - `fg!`: compute cost function and gradient in place
 - `γ`: grid
 """
-function filterdata(u₀,Alu,d₀,y,W⁻,fg!,γ)
+function steadyclimatology(u₀,Alu,d₀,y,W⁻,fg!,γ)
 
     # a first guess: observed surface boundary conditions are perfect.
     # set surface boundary condition to the observations.
@@ -1104,9 +1338,9 @@ end
 - `W⁻`: appropriate weighting (inverse covariance) matrix for these observations,
 - `θtrue`: real observations, 3D field
 """
-function sample_observations(TMIversion,variable)
-    
-    A, Alu, γ, inputfile = config(TMIversion)
+function sample_observations(TMIversion,variable,γ)
+
+    inputfile = datadir("TMI_"*TMIversion*".nc")
 
     # take synthetic observations
     # get observational uncertainty
@@ -1142,9 +1376,9 @@ end
 - `locs`: 3-tuples of locations for observations
 - `wis`: weighted indices for interpolation to locs sites
 """
-function sample_observations(TMIversion,variable,N)
-    
-    A, Alu, γ, inputfile = config(TMIversion)
+function sample_observations(TMIversion,variable,γ,N)
+
+    inputfile = datadir("TMI_"*TMIversion*".nc")
 
     # take synthetic observations
     # get observational uncertainty
@@ -1563,6 +1797,11 @@ end
 function iswet(loc,γ)
     # two approaches
     # approach 1
+
+    # 1 = very strict
+    # 0 = all points
+    wetness = 0.5
+    
     wis = interpindex(loc,γ)
 
     # handle wraparound
@@ -1575,7 +1814,364 @@ function iswet(loc,γ)
     # will be greater than 0.
     # this criterion only requires on land point nearby,
     # where nearby is one of the 8 corners of the cube that contains loc
-    return wetwrap[wis...] > 0.9
+    return wetwrap[wis...] > wetness
+end
+
+""" 
+Save TMI configuration to NetCDF format for non-proprietary access
+"""
+function config2nc(TMIversion,A,γ,L,B)
+
+    # make new netcdf file.
+    filenetcdf = datadir("TMI_"*TMIversion*".nc")
+    isfile(filenetcdf) && rm(filenetcdf)
+
+    grid2nc(TMIversion,γ)
+    
+    matfields2nc(TMIversion,γ)
+
+    !isnothing(A) && watermassmatrix2nc(TMIversion,A)
+
+    !isnothing(L) && circulationmatrix2nc(TMIversion,L,γ)
+
+    !isnothing(B) && boundarymatrix2nc(TMIversion,B)
+
+    #= is this part of the config? Or should it go to
+     a separate output? It is similar to the output fields above. Probably should be considered part of the config. =#
+    regions2nc(TMIversion,γ)
+
+end
+
+"""
+Save grid dictionaries of attributes for writing to NetCDF file
+"""
+function griddicts(γ)
+    # update names and types in dictionary
+    
+    TMIgrids = Dict("lon" => γ.lon,
+                    "lat" => γ.lat,
+                    "depth" => γ.depth)
+    
+    TMIgridsatts = Dict("lon" => Dict("longname" => "Longitude", "units" => "°E"),
+                        "lat" => Dict("longname" => "Latitude", "units" => "°N"),
+                        "depth" => Dict("longname" => "depth", "units" => "m"))
+
+    return TMIgrids, TMIgridsatts
+
+end
+
+"""
+Read 3D fields from mat file and save to NetCDF file.
+"""
+function matfields2nc(TMIversion,γ)
+
+    filenetcdf = datadir("TMI_"*TMIversion*".nc")
+    filemat = datadir("TMI_"*TMIversion*".mat")
+    vars = matread(filemat)
+
+    TMIgrids, TMIgridsatts = griddicts(γ)
+
+    T = eltype(γ.lon) # does the eltype of longitude have to equal the tracer eltype?
+    #T =  Float64
+
+    varlist = Dict("dP" => "qPO₄",
+                   "Tobs" => "θ",
+                   "Terr" => "σθ",
+                   "Sobs" => "Sp",
+                   "Serr" => "σSp",
+                   "O18obs" => "δ¹⁸Ow",
+                   "O18err" => "σδ¹⁸Ow",
+                   "Pobs" => "PO₄",
+                   "Perr" => "σPO₄",
+                   "Nobs" => "NO₃",
+                   "Nerr" => "σNO₃",
+                   "Oobs" =>  "O₂",
+                   "Oerr" =>  "σO₂",
+                   "C13obs" =>  "δ¹³C",
+                   "C13err" =>  "σδ¹³C")
+
+    # iterate over all possible variables listed above
+    Izyx = cartesianindex(filemat)
+    TMIfields = Dict{String,Array{T,3}}()
+    for (kk,vv) in varlist
+        haskey(vars,kk) ? push!(TMIfields, vv => tracerinit(vars[kk], Izyx, γ.wet)) : nothing
+    end
+
+    TMIfieldsatts = fieldsatts()
+
+    # iterate in TMIgrids Dictionary to write to NetCDF.
+    for (varname,varvals) in TMIfields
+        
+        nccreate(filenetcdf,varname,"lon",γ.lon,TMIgridsatts["lon"],"lat",γ.lat,TMIgridsatts["lat"],"depth",γ.depth,TMIgridsatts["depth"],atts=TMIfieldsatts[varname])
+        println("write ",varname)
+        ncwrite(varvals,filenetcdf,varname)
+
+    end
+end
+
+"""
+All variable names and attributes.
+Useful for writing NetCDF files.
+"""
+fieldsatts() = 
+    Dict("θ" => Dict("longname" => "potential temperature", "units" => "°C"),
+         "σθ" => Dict("longname" => "1σ standard error in potential temperature", "units" => "°C"),
+         "Sp" => Dict("longname" => "practical salinity", "units" => "PSS-78"),
+         "σSp" => Dict("longname" => "1σ standard error in practical salinity", "units" => "PSS-78"),
+         "δ¹⁸Ow" => Dict("longname" => "oxygen-18 to oxygen-16 ratio in seawater", "units" => "‰ VSMOW"),
+         "σδ¹⁸Ow" => Dict("longname" => "1σ standard error in oxygen-18 to oxygen-16 ratio in seawater", "units" => "‰ VSMOW"),
+         "PO₄" => Dict("longname" => "phosphate", "units" => "μmol/kg"),
+         "σPO₄" => Dict("longname" => "1σ standard error in phosphate", "units" => "μmol/kg"),
+         "qPO₄" => Dict("longname" => "local source of phosphate", "units" => "μmol/kg"),
+         "NO₃" => Dict("longname" => "nitrate", "units" => "μmol/kg"),
+         "σNO₃" => Dict("longname" => "1σ standard error in nitrate", "units" => "μmol/kg"),
+         "O₂" => Dict("longname" => "dissolved oxygen", "units" => "μmol/kg"),
+         "σO₂" => Dict("longname" => "1σ standard error in dissolved oxygen", "units" => "μmol/kg"),
+         "δ¹³C" => Dict("longname" => "carbon-13 to carbon-12 ratio in DIC", "units" => "‰ PDB"),
+         "σδ¹³C" => Dict("longname" => "1σ standard error fin carbon-13 to carbon-12 ratio in DIC", "units" => "‰ PDB"),
+         "F₀" => Dict("longname" => "normalized mass flux out of gridcell", "units" => "(kg seawater/s)/(kg gridcell)"))
+
+
+"""
+Read vectors from mat file, translate to 3D,
+ and save surface field to NetCDF file.
+"""
+function regions2nc(TMIversion,γ)
+
+    filenetcdf = datadir("TMI_"*TMIversion*".nc")
+    filemat = datadir("TMI_"*TMIversion*".mat")
+
+    # region names
+    # didn't figure out how to use an ordered dict, instead use a tuple
+    list = ("GLOBAL","ANT","SUBANT",
+            "NATL","NPAC","TROP","ARC",
+            "MED","ROSS","WED","LAB","GIN",
+            "ADEL","SUBANTATL","SUBANTPAC","SUBANTIND",
+            "TROPATL","TROPPAC","TROPIND")
+
+    regionname = Dict("GLOBAL" => "globally uniform",
+                      "ANT" => "Antarctic",
+                      "SUBANT" => "Subantarctic",
+                      "NATL" => "North Atlantic",
+                      "NPAC" => "North Pacific",
+                      "TROP" => "tropical and subtropical",
+                      "ARC" => "Arctic",
+                      "MED" => "Mediterranean",
+                      "ROSS" => "Ross Sea sector",
+                      "WED" => "Weddell Sea sector",
+                      "LAB" => "Labrador and Irminger Seas",
+                      "GIN" => "Greenland-Iceland-Norwegian Seas",
+                      "ADEL" => "Adélie Land sector",
+                      "SUBANTATL" => "Atlantic-sector Subantarctic",
+                      "SUBANTPAC" => "Pacific-sector Subantarctic",
+                      "SUBANTIND" => "Indian-sector Subantarctic",
+                      "TROPATL" => "tropical and subtropical Atlantic",
+                      "TROPPAC" => "tropical and subtropical Pacific",
+                      "TROPIND" => "tropical and subtropical Indian")
+    
+    matobj = matopen(filemat)
+    if haskey(matobj,"d_all")
+        d_all = read(matobj,"d_all")
+        close(matobj)
+    else
+        return
+    end
+
+    # a kludge for now
+    T = eltype(γ.lon)
+    
+    # iterate over all regions in d_all
+    Izyx = cartesianindex(filemat)
+    regions = Dict{String,Array{T,2}}()
+    regionatts = Dict{String,Dict{String,String}}()
+    
+    for rr = 1:size(d_all,2)
+        # 3D fields in zyx vector format
+        # are changed to 3D xyz format
+        d = tracerinit(d_all[:,rr],Izyx,γ.wet)
+
+        # just save the surface 2D field
+        push!(regions, list[rr] => d[:,:,1])
+        
+        push!(regionatts, list[rr] =>
+         Dict("longname" => regionname[list[rr]]*" surface region", "units" => "[]"))
+    end
+
+    TMIgrids, TMIgridsatts = griddicts(γ)
+
+    # iterate in regions Dictionary to write to NetCDF.
+    for (varname,varvals) in regions
+        dvarname = "d_"*varname
+        nccreate(filenetcdf,dvarname,"lon",γ.lon,TMIgridsatts["lon"],"lat",γ.lat,TMIgridsatts["lat"],atts=regionatts[varname])
+        println("write ",dvarname)
+        ncwrite(varvals,filenetcdf,dvarname)
+
+    end
+end
+
+
+function watermassmatrix2nc(TMIversion,A)
+
+    filenetcdf = datadir("TMI_"*TMIversion*".nc")
+    i, j, m = findnz(A)
+    nelements = length(i)
+
+    # add the circulation matrix: problem can't store sparse matrix.
+    varname= "m"
+    elementatts = Dict("longname" => "TMI sparse matrix element number")
+    matts =  Dict("longname" => "TMI water-mass fraction (sparse matrix values)", "units" => "(kg source)/(kg total)")
+    nccreate(filenetcdf,varname,"A_element",1:nelements,elementatts,atts=matts)
+    println("write ",varname)
+    ncwrite(m, filenetcdf,varname)
+    
+     varname= "Arow"
+     destatts = Dict("longname" => "gridcell number of destination (row value)")
+     nccreate(filenetcdf,varname,"A_element",1:nelements,elementatts,atts=destatts)
+    println("write ",varname)
+    ncwrite(i, filenetcdf,varname)
+
+    varname= "Acol"
+    sourceatts = Dict("longname" => "gridcell number of source (column value)")
+    nccreate(filenetcdf,varname,"A_element",1:nelements,elementatts,atts=sourceatts)
+    println("write ",varname)
+    ncwrite(j, filenetcdf,varname)
+
+end
+
+"""
+Save circulation matrix `L` to NetCDF file.
+"""
+function circulationmatrix2nc(TMIversion,L,γ)
+
+    T = eltype(L)
+    fullmatrix = false # more efficient to just save F₀, then modify A to get L 
+    filenetcdf = datadir("TMI_"*TMIversion*".nc")
+    if !fullmatrix
+        F₀ = tracerinit(γ.wet,T)
+        for nd ∈ eachindex(γ.I)
+            # normalized mass flux out of gridcell is found on diagonal
+            F₀[γ.I[nd]] = -L[nd,nd]
+        end
+
+        TMIgrids, TMIgridsatts = griddicts(γ)
+        TMIfieldsatts = fieldsatts()
+        varname = "F₀"
+        nccreate(filenetcdf,varname,"lon",γ.lon,TMIgridsatts["lon"],"lat",γ.lat,TMIgridsatts["lat"],"depth",γ.depth,TMIgridsatts["depth"],atts=TMIfieldsatts[varname])
+        println("write ",varname)
+        ncwrite(F₀,filenetcdf,varname)
+
+    else # fullmatrix = true, makes a bigger nc file
+
+        i, j, F = findnz(L)
+        nelements = length(i)
+
+        # add the circulation matrix: problem can't store sparse matrix.
+        varname = "F"
+        elementatts = Dict("longname" => "TMI sparse matrix element number")
+        Fatts =  Dict("longname" => "normalized mass flux (sparse matrix values)","units" => "(kg seawater/s)/(kg gridcell)")
+        nccreate(filenetcdf,varname,"L_element",1:nelements,elementatts,atts=Fatts)
+        println("write ",varname)
+        ncwrite(F, filenetcdf,varname)
+        
+        varname= "Lrow"
+        destatts = Dict("longname" => "gridcell number of destination (row value)")
+        nccreate(filenetcdf,varname,"L_element",1:nelements,elementatts,atts=destatts)
+        println("write ",varname)
+        ncwrite(i, filenetcdf,varname)
+
+        varname= "Lcol"
+        sourceatts = Dict("longname" => "gridcell number of source (column value)")
+        nccreate(filenetcdf,varname,"L_element",1:nelements,elementatts,atts=sourceatts)
+        println("write ",varname)
+        ncwrite(j, filenetcdf,varname)
+    end
+    
+end
+
+"""
+Save boundary matrix for transient model to NetCDF file
+"""
+function boundarymatrix2nc(TMIversion,B)
+
+    filenetcdf = datadir("TMI_"*TMIversion*".nc")
+    i, j, b = findnz(B)
+    nelements = length(i)
+
+    # add the circulation matrix: problem can't store sparse matrix.
+    varname= "b"
+    elementatts = Dict("longname" => "TMI boundary matrix element number")
+    matts =  Dict("longname" => "Boundary matrix values", "units" => "[]")
+    nccreate(filenetcdf,varname,"B_element",1:nelements,elementatts,atts=matts)
+    println("write ",varname)
+    ncwrite(b, filenetcdf,varname)
+
+    varname= "Brow"
+    destatts = Dict("longname" => "gridcell number of 3D field")
+    nccreate(filenetcdf,varname,"B_element",1:nelements,elementatts,atts=destatts)
+    println("write ",varname)
+    ncwrite(i, filenetcdf,varname)
+
+    varname= "Bcol"
+    sourceatts = Dict("longname" => "gridcell number of surface boundary condition")
+    nccreate(filenetcdf,varname,"B_element",1:nelements,elementatts,atts=sourceatts)
+    println("write ",varname)
+    ncwrite(j, filenetcdf,varname)
+
+end
+
+"""
+Put grid properties (Cartesian index) into NetCDF file
+"""
+function grid2nc(TMIversion,γ)
+
+    filenetcdf = datadir("TMI_"*TMIversion*".nc")
+
+    linearindexatts = Dict("longname" => "linear index")
+    nfld = length(γ.I)
+
+    iatts =  Dict("longname" => "Cartesian index for x-direction", "units" => "[]")
+    jatts =  Dict("longname" => "Cartesian index for y-direction", "units" => "[]")
+    katts =  Dict("longname" => "Cartesian index for z-direction", "units" => "[]")
+
+    varname = "i"
+    nccreate(filenetcdf,varname,"linearindex",1:nfld,linearindexatts,atts=iatts)
+    ncwrite(lonindex(γ.I),filenetcdf,varname)
+    
+    varname = "j"
+    nccreate(filenetcdf,varname,"linearindex",1:nfld,linearindexatts,atts=jatts)
+    ncwrite(latindex(γ.I),filenetcdf,varname)
+    
+    varname = "k"
+    nccreate(filenetcdf,varname,"linearindex",1:nfld,linearindexatts,atts=katts)
+    ncwrite(depthindex(γ.I),filenetcdf,varname)
+    
+end
+
+"""
+Read an oceanographically-relevant surface region from NetCDF file. (Also could be read from mat file.)
+"""
+function surfaceregion(TMIversion,region,γ)
+
+    file = datadir("TMI_"*TMIversion*".nc")
+    T = eltype(γ.lon)
+    tracername = "d_"*region
+    dsfc = ncread(file,tracername)
+    println("sumdsfc",sum(filter(!isnan,dsfc)))
+    println("maxdsfc",maximum(filter(!isnan,dsfc)))
+    # expand dsfc to cover 3D.
+    # will not use control2state, because the control
+    # make include non-surface regions in future.
+
+    # preallocate
+    d = Array{T}(undef,size(γ.wet))
+
+    # set ocean to zero, land to NaN
+    # consider whether land should be nothing or missing
+    d[γ.wet] .= zero(T)
+    d[.!γ.wet] .= zero(T)/zero(T)
+    d[:,:,1] = dsfc
+    println(count(!isnan,d))
+    return d
 end
 
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 1f34bda4..9e0ecf67 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,8 +2,12 @@ using Revise, TMI, Test
 
 @testset "TMI.jl" begin
 
-    TMIversion = "TMI_2010_2012_4x4x33"
-    A, Alu, γ, TMIfile = config(TMIversion)
+    TMIversion = "modern_90x45x33_GH10_GH12"
+    #TMIversion = "modern_180x90x33_GH10_GH12"
+    #TMIversion = "modern_90x45x33_unpub12"
+    #TMIversion = "modern_90x45x33_G14"
+    
+    A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)
     
     ############################
     ## trackpathways
@@ -13,12 +17,12 @@ using Revise, TMI, Test
         @test minimum(sum(A,dims=2))> -1e-14
         
         #- define the surface patch by the bounding latitude and longitude.
-        latbox = [50,60]; # 50 N -> 60 N, for example.
+        latbox = [50. , 60.]; # 50 N -> 60 N, for example.
 
         # mutable due to wraparound: don't use an immutable tuple
-        lonbox = [-50,0]; # 50 W -> prime meridian
+        lonbox = [-50. , 0.]; # 50 W -> prime meridian
 
-        c,γ = trackpathways(TMIversion,latbox,lonbox)
+        c = trackpathways(Alu,latbox,lonbox,γ)
 
         @test maximum(c[γ.wet]) ≤ 1.0
         @test minimum(c[γ.wet]) ≥ 0.0
@@ -28,7 +32,7 @@ using Revise, TMI, Test
     ## example: regeneration
     @testset "regeneration" begin
 
-        PO₄ᴿ, γ = regeneratedphosphate(TMIversion)
+        PO₄ᴿ = regeneratedphosphate(TMIversion,Alu,γ)
         @test maximum(PO₄ᴿ[γ.wet]) < 10
         @test minimum(PO₄ᴿ[γ.wet]) ≥ 0
     end
@@ -46,7 +50,7 @@ using Revise, TMI, Test
         dVdd = tracerinit(γ.wet); # pre-allocate c
         dVdd[γ.wet] = Alu'\v[γ.wet]
 
-        volume = volumefilled(TMIversion)
+        volume = volumefilled(TMIversion,Alu,γ)
 
         # volumefill positive at surface?
         @test sum(volume .< 0) == 0
@@ -60,11 +64,7 @@ using Revise, TMI, Test
         # choose linear or nearest neighbor interpolation
         linearinterp = true #true
         
-        # would be better to randomize location.
-        # xlon = 125.26; # deg E.
-        # xlat = -6.38;  # deg N.
-        # xdepth = 2578.2;  # meters.
-        # loc = (xlon,xlat,xdepth)
+        # randomized location.
         loc = wetlocation(γ)
 
         # choose linear or nearest neighbor interpolation
@@ -80,7 +80,7 @@ using Revise, TMI, Test
 
         @test isapprox(sum(filter(!isnan,δ)),1.0) 
 
-        origin, γ = surfaceorigin(TMIversion,loc)
+        origin = surfaceorigin(loc, Alu, γ)
         
         @test isapprox(sum(filter(!isnan,origin)),1)
         @test isapprox(sum(filter(!isnan,origin)),1)
@@ -89,15 +89,15 @@ using Revise, TMI, Test
     end
     
     #########################################
-    ## filterdata
-    @testset "filterdata" begin
+    ## formerly filterdata
+    @testset "steadyclimatology" begin
 
         # first guess of change to surface boundary conditions
         # ocean values are 0
         u₀ = zeros(Float64,sum(γ.wet[:,:,1]))
 
         # take synthetic, noisy observations
-        y, W⁻, ctrue = sample_observations(TMIversion,"θ")
+        y, W⁻, ctrue = sample_observations(TMIversion,"θ",γ)
 
         # a first guess: observed surface boundary conditions are perfect.
         # set surface boundary condition to the observations.
@@ -113,7 +113,7 @@ using Revise, TMI, Test
         #fg!(J̃₀,gJ₀,u₀)
         # filter the data with an Optim.jl method
 
-        out = filterdata(u₀,Alu,y,d₀,W⁻,fg!,γ)
+        out = steadyclimatology(u₀,Alu,y,d₀,W⁻,fg!,γ)
 
         # check with forward differences
         ϵ = 1e-5
@@ -149,7 +149,7 @@ using Revise, TMI, Test
         u₀ = zeros(Float64,sum(γ.wet[:,:,1]))
 
         # take synthetic, noisy observations
-        y, W⁻, ctrue, locs, wis = sample_observations(TMIversion,"θ",N)
+        y, W⁻, ctrue, locs, wis = sample_observations(TMIversion,"θ",γ,N)
 
         # does this help optimization stay stable?
         #W⁻ *= 1.0/100.0
@@ -157,7 +157,8 @@ using Revise, TMI, Test
         # make a silly first guess for surface
         d₀ = tracerinit(γ.wet)
         [d₀[γ.I[ii]] = 15.0 for ii ∈ eachindex(γ.I) if γ.I[ii][3] == 1]
-                
+
+        # permit surface deviations on order of 5°C
         Q⁻ = 1.0/(5.0^2)
         #Q⁻ = 10.0