TMI

struct Field

This structure permits the grid to be 
automatically passed to functions with
the tracer field.

This structure assumes the Tracer type to be 
three-dimensional.

tracer::AbstractArray{T,N}
γ::Grid{A,N}
name::Symbol
longname::String
units::String

struct Grid

TMI grid with accounting for wet/dry points

Fields

• axes::NTuple{N,Vector{A}}: labels for axis such as lon, lat, depth with element type A
• wet::BitArray{N}: mask for ocean points
• interior::BitArray{N}: mask for interior ocean points
• wrap::NTuple{N,Bool}: does the domain wraparound in each dim?
• Δ::Vector{CartesianIndex{N}}: defines computational stencil relative to central cell

function Grid(TMIfile)

Construct the Grid given a file name

Arguments

• TMIfile::String: NetCDF file name for TMI version

Output

• γ::Grid: TMI grid struct

function Grid(foreign_file, maskname, lonname, latname, depthname)

Construct the Grid from a non-TMI file given the names of relevant fields.

Assumes that an ocean mask is available.
Assumes an input NetCDF file.
Assumes everything below the top layer is part of the interior. 
Tested for Float32 fields (should work for other types).

Arguments

• foreign_file::String
• maskname::String
• lonname::String
• latname::String
• depthname::String

Output

• γ::Grid: TMI grid struct

struct MassFraction

store mass fractions in a Field-like array

• fraction::Array{T,3}
• γ::Grid
• name::Symbol
• longname::String
• units::String
• position::CartesianIndex{3}

struct Source

This structure describes Sources, which are
similar to Fields, but they may be
1) non-negative
2) have only interior mask

`function \(A,d::Field)::Field`
Define left division for Fields
Need two slashes to prevent invalid escape

Base.length(c::Union{Field,Source,BoundaryCondition}) = length(c.tracer[wet(c)])

Extend `length` to give the number of wet (i.e., ocean) gridcells.

function oneunit, help for gridded Interpolations

function oneunit, help for gridded Interpolations

function ones(γ::Grid,name=:none,longname="unknown",units="unknown")::Field

  initialize tracer field of ones on TMI grid
  using a Field struct and constructor

Arguments

• γ::TMI.Grid

Output

• d::Field, 3d tracer field with NaN on dry points

function ones(dim::Int64,dimval::Int64,γ::Grid)::BoundaryCondition

   Initialize boundary condition with ones

Base.propertynames(γ::Grid) = (I,R,fieldnames(typeof(x))...)

Do not store Cartesian and linear indices.
Compute them on demand.

Specialize Base.sum(c::Field)

so that it doesn't use the slow iteration method

function vec(u)

Turn a collection of controls into a vector
for use with Optim.jl. 
An in-place version of this function would be handy.

zero(c::Field) = zeros(c.γ)

function zeros(wet,ltype=Float64)
initialize tracer field on TMI grid
This version will give an array

Arguments

• wet::BitArray mask of ocean points
• ltype:: optional type argument, default=Float64

Output

• d:: 3d tracer field with NaN on dry points

function zeros(γ::Grid,name=:none,longname="unknown",units="unknown")::Field

  initialize tracer field on TMI grid
  using a Field struct and constructor

Arguments

• γ::TMI.Grid

Output

• d::Field, 3d tracer field with NaN on dry points

function zeros(dim::Int64,dimval::Int64,γ::Grid,name::Symbol,longname::String,units::String)::BoundaryCondition

   Initialize boundary condition with zeroes

Arguments

• dim:
• dimval
• γ::Grid
• name::Symbol
• longname::String
• units::String

Output

• b::BoundaryCondition

`function *(c::Field,d::Field)::Field`
Field by field multiplication is element-by-element.

function Statistics.mean(c::Field)

Take the volume-weighted mean of a `Field`

`function +(c::BoundaryCondition,d::BoundaryCondition)::BoundaryCondition`
Define addition for Fields

function adjustboundarycondition!(b::Union{BoundaryCondition,NamedTuple},u::Union{BoundaryCondition,NamedTuple})

adjust all boundary conditions b that are described in u

warning: if u doesn't contain any boundary condition adjustments,
nothing will change.

function adjustboundarycondition!(b::Union{BoundaryCondition,NamedTuple},u::Union{BoundaryCondition,NamedTuple})

adjust all boundary conditions b that are described in u

function axislabels(file)
Read and assemble the grid properties.

Arguments

• file: TMI NetCDF file name

Output

• grid: TMI grid coordinates

function boundaryconditionatts(dim::Int64,dimval::Int64,γ::Grid)

   Help initialize boundary condition by getting some attributes

    function boundarymatrix(file,γ)
Read and assemble the boundary matrix from MATLAB.
Transfer to updated x,y,z version

Arguments

• file: TMI MATLAB file name
• γ: TMI grid

Output

• B: boundary condition matrix

Save boundary matrix for transient model to NetCDF 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

• I: TMI Cartesian index for wet points

function cartesianindex(wet)
Read and assemble the grid coordinates
according to a 3D tracer in x,y,z order

Arguments

• wet: BitArray logical mask for wet points

Output

• I: 3D Cartesian indices

Horizontal area of grid cell

function cellvolume(γ)::Field

Volume of each grid cell.

function checkgrid!(c,wet)

perform a check of file compatibility
 with grid

function circulationmatrix(file,A,γ)
Read and assemble the circulation matrix from the efficient storage of A and F₀ variables.

Arguments

• file: TMI MATLAB file name
• A: TMI water-mass matrix
• γ: TMI grid

Output

• L: circulation matrix in xyz format

function circulationmatrix(file,γ)
Read and assemble the circulation matrix from MATLAB.
Transfer to updated x,y,z version

Arguments

• file: TMI MATLAB file name
• γ: TMI grid

Output

• L: circulation matrix in xyz format

Save circulation matrix L to NetCDF file.

Save TMI configuration to NetCDF format for non-proprietary access

Configure TMI environment from original MATLAB output

function config_from_nc(TMIversion; compute_lu = true)
Configure TMI environment from NetCDF input file.

Arguments

• TMIversion: TMI version for water-mass/circulation model

Output

• A: TMI steady-state water-mass matrix
• Alu: LU decomposition of A
• γ: TMI grid properties
• TMIfile: TMI file name

function control2state(tracer2D,γ)
turn 2D surface field into 3D field with zeroes below surface

Arguments

• tracer2D:: 2D surface tracer field
• wet::BitArray mask of ocean points

Output

• tracer3D:: 3d tracer field with NaN on dry points

function costfunction_gridded_model!(J,guvec,convec::Vector{T},non_zero_indices,u₀::Union{BoundaryCondition{T},NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}},c,y::Field{T},Wⁱ::Diagonal{T, Vector{T}},Qⁱ::Diagonal{T, Vector{T}},γ::Grid) where {N1, N2, T <: Real}

function costfunction_gridded_model(convec::Vector{T},non_zero_indices,y::Field{T},u,A0,c,q,Wⁱ::Diagonal{T, Vector{T}},Qⁱ::Diagonal{T, Vector{T}},γ::Grid) where T <: Real

squared model-data misfit for gridded data
controls are a vector input for Optim.jl

Arguments

• convec: concatenated control vecotr incuding u and f
• J: cost function of sum of squared misfits
• gJ: derivative of cost function wrt to controls
• u: tracer controls, field format
• non_zero_indices: Non-zero indices for reconstruction of water-mass matrix A
• c: tracer concentrations from GCM
• Wⁱ: inverse of W weighting matrix for observations
• Qⁱ: inverse of Q weighting matrix for tracer conservation
• γ: grid

function costfunction_point_obs!(J,guvec,uvec::Vector{T},Alu,b₀::BoundaryCondition{T},u₀::BoundaryCondition{T},y::Vector{T},Wⁱ::Diagonal{T, Vector{T}},wis,locs,Q⁻,γ::Grid) where T <: Real

squared model-data misfit for pointwise data
controls are a vector input for Optim.jl
Issue #1: couldn't figure out how to nest with costfunction_obs!

Arguments

• J: cost function of sum of squared misfits
• guvec: derivative of cost function wrt to controls
• uvec: controls, vector format
• Alu: LU decomposition of water-mass matrix
• b: boundary condition
• y: pointwise observations
• Wⁱ: inverse of W weighting matrix for observations
• wis: weights for interpolation (data sampling, E)
• locs: data locations (lon,lat,depth)
• Q⁻: weights for control vector
• γ: grid

function costfunction_point_obs(uvec::Vector{T},Alu,b₀::BoundaryCondition{T},u₀::BoundaryCondition{T},y::Vector{T},Wⁱ::Diagonal{T, Vector{T}},wis,locs,Q⁻,γ::Grid;q=nothing,r=1.0) where T <: Real

Squared model-data misfit for pointwise data.
Controls are a vector input for Optim.jl.
Core numerics handled by `costfunction_point_obs`.

Arguments

• uvec: controls, vector format
• Alu: LU decomposition of water-mass matrix
• b: boundary condition
• y: pointwise observations
• Wⁱ: inverse of W weighting matrix for observations
• wis: weights for interpolation (data sampling, E)
• locs: data locations (lon,lat,depth)
• Q⁻: weights for control vector
• γ: grid

Optional

• q::Field: interior source
• r::Number: scalar factor for source

Output

• J: cost function of sum of squared misfits
• gJ: derivative of cost function wrt to controls

function depthindex(I) 
Get the k-index (depth level) from the Cartesian index

function distancematrix(γ;surface = true)

Matrix with size of surface points squared

Each entry gives distance in km between surface points
Gives only horizontal distance.

function download_matfile(TMIversion::String)

Download MATLAB file for given TMI version

Arguments

• TMIversion::String

Output

• TMIfile::String: TMI file name

Side-effect

• download TMI input file if necessary

function download_ncfile(TMIversion::String)

Download NetCDF file for given TMI version

Arguments

• TMIversion::String

Output

• TMIfile::String: TMI file name

Side-effect

• download TMI input file if necessary

function download_regionfile(TMIversion::String)

Return file name of regional masks.

Also download file from Google Drive, if not already done.

File name refers to the 2D size of domain. It is the same for modern and LGM configs and only depends on number of points in Nx and Ny directions.

function eastindex(I) 
Get the vector index on the northern open boundary

All variable names and attributes. Useful for writing NetCDF files.

function gadjustboundarycondition!(b::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real

adjust the (one) boundary condition 
Just copy the variable.
Keep this function so that calling functions can look alike.
Could probably combine with lower function, use Union type

function gaussiandistancematrix(γ,σ,L)

uses distance matrix plus a lengthscale `L` (km)
and a size of the diagonal `σ`

returns values with inverse gaussian weighting

function getboundarycondition(field::Field,dim,dimval)::BoundaryCondition

Get boundary condition by extracting from Field (i.e., 3D tracer)

Arguments

• field::Field: 3D tracer field with metadata and grid
• dim: dimension number (1,2,3) that the boundary plane has constant value
• dimval: index number in dimension dim that defines boundary plane

Output

• b::BoundaryCondition: boundary condition on a plane with metadata and grid

Get boundary condition by extracting from N-dimensional tracer and returning (N-1)-dimensional array

function gobserve(gy::Vector{T},c::Field{T},wis,γ) where T <: Real

ADJOINT Take a observation at location given by weights wis
Arguments not symmetric with `observe` due to splat operator

Put grid properties (Cartesian index) into NetCDF file

Save grid dictionaries of attributes for writing to NetCDF file

function gridsize(TMIversion)

Parse the TMIversion string for the grid size

Will break if the prefix contains underscores or 'x'

function gsetboundarycondition(gd::Field{T},b::BoundaryCondition{T}) where T<: Real

ADJOINT: apply boundary condition to the equation constraints

Arguments

• d::Field, equation constraints (i.e., right hand side)
• b::BoundaryCondition

function gsetboundarycondition(gd::Field{T},b::BoundaryCondition{T}) where T<: Real

ADJOINT: apply boundary condition to the equation constraints

Arguments

• d::Field, equation constraints (i.e., right hand side)
• b::BoundaryCondition

function gsetsource!(gq::Field{T},gd::Field{T},r=1.0)

Adjoint to `setsource!`

function gsteadyinversion(gc,Alu,b;q=nothing,r=1.0)

ADJOINT invert for a steady-state tracer distribution

Arguments

• Alu: LU decomposition of water-mass matrix
• b: BoundaryCondition
• γ::Grid

Optional Arguments

• q: interior sources/sinks of phosphate
• r: stochiometric ratio of tracer:phosphate

Output

• c::Field, steady-state tracer distribution

function horizontaldistance(loc,γ)

Arguments

• loc: 3-tuple of lon,lat,depth location
• γ: TMI.grid

Output

• hordist: horizontal distance to nearest tracer grid points

function interiormask(A,wet,nx,ny,nz)

function interpindex(loc,γ) Weights for linear interpolation. The derivative of linear interpolation is needed in sensitivity studies. ReverseDiff.jl could find this quantity automatically. Instead we dig into the Interpolations.jl package to find the weights that are effectively the partial derivatives of the function.

Arguments

• c: a temporary tracer field, would be nice to make it unnecessary
• loc: (lon,lat,depth) tuple of a location of interest
• γ: TMI grid

Output

• δ: weights on a 3D tracer field grid

function interpweights(loc,γ) Weights for linear interpolation. The derivative of linear interpolation is needed in sensitivity studies. ReverseDiff.jl could find this quantity automatically. Instead we dig into the Interpolations.jl package to find the weights that are effectively the partial derivatives of the function.

Arguments

• loc: (lon,lat,depth) tuple of a location of interest
• γ: TMI grid

Output

• δ: weights on a 3D tracer field grid

function latindex(I) 
Get the j-index (latitude index) from the Cartesian index

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 local_solve(c::NamedTuple, w::NamedTuple)

Given tracers, solve for mass fractions using a repeated local algorithm.

local_solve! is the workhorse for this algorithm and specifies a default inverse tapering parameter of α=100_000.

Arguments

• c::NamedTuple: Field tracers (observations or model output)
• w::NamedTuple: scale size of tracers (for weighting observational fits)

Output

• m̃::NamedTuple: MassFraction in a NamedTuple collection

function local_watermass_matrix(c::NamedTuple, m::NamedTuple, I::CartesianIndex, neighbors::Field) Find local water-mass matrix with singularity checker (true if one neighbor only has a single connection to the rest of the ocean)

Arguments

• c::NamedTuple: input tracers
• m::NamedTuple: mass fractions for grid stencil
• I::CartesianIndex: local "location"
• neighbors::Field: integer number of neighbors

Output

• A::Matrix: local water-mass matrix
• single_connection::Bool: true if flagged for singularity warning

function lonindex(I) 
Get the i-index (lon index) from the Cartesian index

function massfractions(c::NamedTuple, w::NamedTuple; alg = :local)

Create NamedTuple of mass fractions from observations c

Doesn't produce a MassFraction struct and thus is named in lower case.

Arguments

• c::NamedTuple: input observations
• w::NamedTuple: scale size of observations (used if obs not fit exactly)
• alg=:local: default algorithm is :local

Output

• m::NamedTuple: collection of mass fractions

function mat2ncfield

Rename MATLAB variables to NetCDF variables

Read 3D fields from mat file and save to NetCDF file.

Read 3D fields from mat file and save to NetCDF file.

    function matrix_zyx2xyz(TMIfile,Azyx,γ)

Transfer zyx format water-mass matrix A to xyz format

Arguments

• Azyx: water-mass matrix in zyx format
• γ: TMI grid

Output

• Axyz: water-mass matrix in xyz format

Read 3D source field from mat file and save to NetCDF file.

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 meanage(TMIversion,Alu,γ)
Mean or ideal age

Arguments

• TMIversion: version of TMI water-mass/circulation model
• Alu: LU decomposition of water-mass matrix A
• γ: TMI grid

Output

• a: mean age [yr]

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 nearestneighbor(loc,γ)
return the Cartesian index and linear index 
of the nearest N neighbors

Arguments

• loc: 3-tuple of lon,lat,depth location
• γ: TMI.grid

Output

• Inn: Cartesian indices of nearest neighbor

#- Rnn: linear indices of nearest neighbor, Removed from code

function nearestneighbormask
Make a 3D tracer field that is 1 at location 
of nearest neighbor, 0 elsewhere

Arguments

• loc: location in a 3-tuple (lon,lat,depth)
• γ: TMI.grid

Output

• δ: nearest neighbor mask 3D field

function neighbor_indices(n::Integer)

Direction (step) of neighbors away from a central point. Choose n = 6 (default) or n=26.

Argument

• n=6: max number of neighbors

Output

• In::Vector{CartesianIndex}: indices of neighbors

function neighbors(γ::Grid; longname = "number of neighbors")

How many neighbors does each grid cell have? Conceptually, it only depends on the grid, but this algorithm is slower than the one that takes mass fractions as input.

Arguments

• γ::TMI.Grid

Output

• n::Field: integer number of neighbors

function neighbors(m::NamedTuple,γ::Grid)

How many neighbors does each grid cell have?

Arguments

• m::NamedTuple: input mass fractions to obtain their stencil (opportunity to simplify)
• γ::TMI.Grid

Output

• n::Field: integer number of neighbors

function northindex(I) 
Get the vector index on the northern open boundary

function observe
Take a observation at location given by weights wis

function observe(c,loc,γ)

Extend the TMI.observe method to use locations rather than weighted interpolations.

Save optimization parameters to NetCDF file)

Future considerations: split into 2 functions

1. read from mat
2. save to nc

preformedcarbon13(TMIversion,Alu,γ) = preformednutrient("δ¹³C",TMIversion,Alu,γ)

preformednitrate(TMIversion,Alu,γ) = preformednutrient("NO₃",TMIversion,Alu,γ)

function preformednutrient(tracer::Union{String,Symbol},TMIversion,Alu,γ)

Preformed (i.e., NO accumulation or remineralization) nutrient

Arguments

• tracer::Union{Symbol,String}: tracer name
• TMIversion: version of TMI water-mass/circulation model
• Alu: LU decomposition of water-mass matrix A
• γ: TMI grid

Output

• c★: preformed tracer

preformedoxygen(TMIversion,Alu,γ) = preformednutrient("O₂",TMIversion,Alu,γ)

preformedphosphate(TMIversion,Alu,γ) = preformednutrient("PO₄",TMIversion,Alu,γ)

function readfield(file,tracername,γ)
Read a tracer field from NetCDF but return it 
as a Field.

Use NCDatasets so that Unicode is correct

Arguments

• file: TMI NetCDF file name
• tracername: name of tracer
• γ::Grid, TMI grid specification

Output

• c::Field

MATLAB version
function readfield(matfile,mattracername,γ::Grid,Izyx) # for MATLAB

read MATLAB field and transfer zyx format to xyz

function readtracer(file,tracername)
Read a tracer field from NetCDF.

Arguments

• file: TMI NetCDF file name
• tracername: name of tracer

Output

• c: 3D tracer field

function regeneratednutrient(TMIversion,Alu,γ)

Regenerated (i.e., accumulated, remineralized) nutrient

Arguments

• tracer::Union{String,Symbol}: tracer name
• TMIversion: version of TMI water-mass/circulation model
• Alu: LU decomposition of water-mass matrix A
• γ: TMI grid
• r: optional stoichiometric ratio relative to PO₄

Output

• PO₄ᴿ: regenerated phosphate

function regions2nc(TMIversion,γ)

Read vectors from mat file, translate to 3D, and save surface field to NetCDF file.

Consider deprecating this function.

function regionurl(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

• file: name of file to look for on Google Drive

Output

• regionurl: location (URL) for download of regional mask file

function section
View latitude-depth slice of field

Arguments

• c::Field, 3D tracer field plus meta data
• lon: longitude of section

Output

• csection: 2d slice of field

function setboundarycondition!(d::Field,b::BoundaryCondition)
apply boundary condition to the equation constraints

Arguments

• d::Field, equation constraints (i.e., right hand side)
• b::BoundaryCondition

function setboundarycondition!(d::Field{T},b::NamedTuple{<:Any, NTuple{N,BoundaryCondition{T}}}) where {N, T <: Real}

set all boundary conditions in a Named Tuple

function setsource!(d::Field,q::Field,r::Number)
apply interior source q to the equation constraints d

Arguments

• d::Field, equation constraints (i.e., right hand side)
• q::Field, interior source
• r::Number, default = 1.0, stoichiometric ratio

function shiftloc(loc)

sometimes loc longitudes are outside of grid due to different conventions
assumption: 360° shift is enough to get back on grid

function southindex(I) 
Get the vector-index on the southern open boundary

function sparsedatamap(u₀::Vector{T},Alu,b::BoundaryCondition{T},u::BoundaryCondition{T},y::Vector{T},W⁻,wis,locs,Q⁻,γ::Grid;iterations=10) where T <: Real

 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₀,
 (c) perturbations to the surface boundary condition u₀
that best fits observations, y,
according to the cost function,
J = (ỹ - y)ᵀ W⁻¹ (ỹ - y)
subject to Aỹ = d₀ + Γ u₀.                 
W⁻ is a (sparse) weighting matrix.
See Supplementary Section 2, Gebbie & Huybers 2011.

Arguments

• u₀:
• Alu:
• b: first guess of boundary conditions and interior sources
• y: observations on 3D grid
• W⁻: weighting matrix best chosen as inverse error covariance matrix
• fg!: compute cost function and gradient in place
• γ: grid

function steadyclimatology_optim(u₀,fg!,iterations) 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₀, (c) perturbations to the surface boundary condition u₀ that best fits observations, y, according to the cost function, J = (ỹ - y)ᵀ W⁻¹ (ỹ - y) subject to Aỹ = d₀ + Γ u₀. W⁻ is a (sparse) weighting matrix. See Supplementary Section 2, Gebbie & Huybers 2011.

Arguments

• u₀:
• fg!: compute cost function and gradient in place
• iterations: number of optimization iterations

function steadyinversion(Alu,b;q=nothing,r=1.0)
invert for a steady-state tracer distribution

Arguments

• Alu: LU decomposition of water-mass matrix
• b: boundary condition, assumed to be surface boundary condition
• γ::Grid

Optional Arguments

• q: interior sources/sinks of phosphate
• r: stochiometric ratio of tracer:phosphate

Output

• c::Field, steady-state tracer distribution

function step_cartesian(I, Δ, γ)

Arguments

• I::CartesianIndex: starting point
• Δ::CartesianIndex: step
• γ::Grid: TMI-defined grid

Output

• Istep::CartesianIndex: new location
• inbounds::Bool: inside the domain bounds?

function surfacecontrol2field!(c,u,γ)
Add surface control vector to existing 3D field

Arguments

• c:: state field, 3d tracer field with NaN on dry points, modified by function
• usfc:: surface control vector
• wet::BitArray mask of ocean points

function surfacecontrol2field!(c,u,γ)
Add surface control vector to tracer vector

Arguments

• c:: state field, 3d tracer field with NaN on dry points, modified by function
• u:: surface control vector
• wet::BitArray mask of ocean points

function surfacecontrol2field(usfc,γ.wet)
turn surface control vector into 3D field with zeroes below surface

Arguments

• usfc:: surface control vector
• wet::BitArray mask of ocean points

Output

• tracer3D:: 3d tracer field with NaN on dry points

function surfaceindex(I) 
Get the vector-index where depth level == 1 and it is ocean.

function surfaceorigin(TMIversion,loc)
 Find the surface origin of water for some interior box 
 This is equivalent to solving a sensitivity problem:
 The mass fraction at a location `loc` of interest is 
`c[loc] = δᵀ c`, where `δ` samples the location of the global mass-fraction variable, c.
Then the sensitivity of `c[loc]` is: d(c[loc])/d(d) = A⁻ᵀ δ.
The derivative is solved using the constraint: Ac = d.
The sensitivity is exactly the mass fraction originating from each source.      
This problem is mathematically similar to determining how the ocean is filled.

Arguments

• 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, log10 of effective depth, in terms of a BoundaryCondition

function surfacepatch
Make a surface boundary condition
with a rectangular patch

Arguments

• lonbox: longitudes of box edges
• latbox: latitudes of box edges
• γ: TMI.grid

Output

• d: vector that describes surface patch

 function surfaceregion(TMIversion::String,region::String,γ::Grid)::BoundaryCondition

Read an oceanographically-relevant surface region from NetCDF file. (Also could be read from mat file.)

Return a BoundaryCondition

Version 1: operates on a 2D Float field

 function surfaceregion(TMIversion::String,region::String,γ::Grid)::BoundaryCondition

Read an oceanographically-relevant surface region from NetCDF file. (Also could be read from mat file.)
Return a BoundaryCondition

function synthetic_observations(TMIversion,variable,locs)
Synthetic observations that are a contaminated version of real observations
This version: observations with random (uniform) spatial sampling

Arguments

• TMIversion::String: version of TMI water-mass/circulation model
• variable::String: variable name to use as template
• N: number of observations

Output

• y: contaminated observations on 3D grid
• W⁻: appropriate weighting (inverse covariance) matrix for these observations,
• ytrue: uncontaminated observations, 3D field
• locs: 3-tuples of locations for observations
• wis: weighted indices for interpolation to locs sites

function synthetic_observations(TMIversion,variable)
Synthetic observations that are a contaminated version of real observations
This version: gridded observations

Arguments

• TMIversion::String: version of TMI water-mass/circulation model
• variable::String: variable name to use as template

Output

• y: contaminated observations on 3D grid
• W⁻: appropriate weighting (inverse covariance) matrix for these observations,
• θtrue: real observations, 3D field

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 trackpathways(TMIversion,latbox,lonbox)
Track the pathways of a user-defined water mass.
 Steps: (a) define the water mass by a rectangular surface patch dyed with passive tracer concentration of         (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
• latbox: min and max latitude of box
• lonbox: min and max longitude of box
• γ: TMI grid

Output

• c: fraction of water from surface source

function unvec!(u,uvec)

Undo the operations by vec(u)
Needs to update u because attributes of 
u need to be known at runtime.

function unvec(u,uvec)

Replace u with new u
Undo the operations by vec(u)
Needs to update u because attributes of 
u need to be known at runtime.

function updatelinearindex(izyx,Izyx,R)
Linear index translated from z,y,x to x,y,z accounting

get Izyx Cartesian index stored from legacy MATLAB code

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

Output

• ixyz: index of interest in x,y,z accounting

function vec2fld
Transfer a vector to a 3D field with accounting for ocean bathymetry

Arguments

• vector: field in vector form (no land points)
• I: cartesian indices of ocean points

Output

• field: field in 3d form including land points (NaN)

function volumefilled(TMIversion)
Find the ocean volume that has originated from each surface box.
 This is equivalent to solving a sensitivity problem:
 The total volume is V = vᵀ c , where v is the volume of each box 
 and c is the fraction of volume from a given source which
 satisfies the equation A c = d.                     
 Next, dV/d(d) = A⁻ᵀ v, and dV/d(d) is exactly the volume originating from each source.

 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: log10 of global ocean volume filled by a surface region, exists at surface, therefore given BoundaryCondition type

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 watermassmatrix(file)
Read and assemble the water-mass matrix.

Arguments

• file: TMI NetCDF or MATLAB file name

Output

• A: water-mass matrix

function watermassmatrix(m::Union{NamedTuple,Vector}, γ::Grid)

Produce water-mass matrix from mass fractions and grid.

Arguments

• m::NamedTuple: collection of MassFractions
• γ::TMI.Grid

Output

• A: sparse water-mass matrix

function westindex(I) 
Get the vector index on the western open boundary

function wetlocation(γ)
Get (lon,lat,depth) tuples of wet locations.
Allow a location to be wet if at least one out of 8 nearby gridpoints is wet.
Certainly "wet" gridpoints could be defined more strictly.

Arguments

• γ: TMI.grid

Output

• loc: lon,lat,depth

function writefield(file,field)

Write a Field to NetCDF.

Use NCDatasets so that Unicode is correct

Arguments

• file: TMI NetCDF file name
• field::Field: a TMI.Field struct

Output

• none

Side-effect

• write to file

function write(file,b)

Write a BoundaryCondition to NetCDF.

Use NCDatasets so that Unicode is correct

Arguments

• file: TMI NetCDF file name
• b::BoundaryCondition: a TMI.BoundaryCondition struct

Output

• none

Side-effect

• write to file

function Γsfc! 
Γsfc! anonymously calls surfacecontrol2field!