TMI

struct BoundaryCondition

a plane defined at `dim=dimval`

Attributes

`tracer::Array{T,2}`: values on plane
`i::Vector{T}`: coordinate values on local x-plane
`j::Vector{T}`: coordinate values on local y-plane
`k::T`: fixed coordinate value on local z-plane that defines the Boundary Condition plane
`dim::Int64`: dimension (1,2, or 3) along which the plane's index is fixed
`dimval::Int64`: plane defined at dim = dimval where dimval is a 1-based index number
`wet::BitArray{2}`: ocean mask for boundary condition

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 +(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 boundarymask(b, γ::Grid)

function boundarymask(γ)

    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 NetCDF.

Arguments

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

Output

• L: circulation matrix in xyz format

Save circulation matrix L to NetCDF file.

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

Save TMI configuration to NetCDF format for non-proprietary access

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_gridded_obs!(J,guvec,uvec::Vector{T},Alu,b₀::Union{BoundaryCondition{T},NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}}},u₀::Union{BoundaryCondition{T},NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}},y::Field{T},Wⁱ::Diagonal{T, Vector{T}},γ::Grid) where {N1, N2, T <: Real}

function costfunction_gridded_obs(uvec::Vector{T},Alu,b::BoundaryCondition{T},y::Field{T},Wⁱ::Diagonal{T, Vector{T}},γ::Grid) where T <: Real

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

Arguments

• J: cost function of sum of squared misfits
• gJ: derivative of cost function wrt to controls
• u: controls, field format
• Alu: LU decomposition of water-mass matrix
• b: boundary conditions
• y: observations on grid
• Wⁱ: inverse of W weighting matrix for observations
• γ: 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 dirichletmatrix(γ, τ)

Dirichlet surface boundary matrix with uniform timescale. Assumes that the Dirichlet boundary condition is zero.

Arguments

• γ: TMI grid
• τ: uniform restoring timescale (years) for all boundary points

Output

• Ldir: circulation matrix in xyz format for boundary points

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_file(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

 function effective_endmember(TMIversion,Alu,field,region,γ)

Effective (i.e., importance-weighted) endmember tracer value calculated according to Gebbie and Huybers 2011.

 function effective_endmember_sums(Alu,field::Field,b::BoundaryCondition,γ::Grid)

Intermediate quantities for computing effective endmembers

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 gadjustboundarycondition(gb::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real

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

geteastboundary(c::Field) = getboundarycondition(c,1,eastindex(c.γ))::BoundaryCondition

getnorthboundary(c::Field) = getboundarycondition(c,2,northindex(c.γ))::BoundaryCondition

getsouthboundary(c::Field) = getboundarycondition(c,2,southindex(c.γ))::BoundaryCondition

getsurfaceboundary(c::Field) = getboundarycondition(c,3,surfaceindex(c.γ))::BoundaryCondition

getwestboundary(c::Field) = getboundarycondition(c,1,westindex(c.γ))::BoundaryCondition

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 iswet(loc, γ, neighbors)
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

• loc: lon,lat,depth
• γ
• neighbors

Output

• Boolean (true for ocean point)

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_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 location_obs(field, locs, γ)

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

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

    function matrix_mod2glacial(Amodern,γmodern,γglacial)

Transfer modern water-mass matrix Amod to glacial water-mass matrix Aglacial

Arguments

• Amodern: modern water-mass matrix
• γmodern: modern TMI grid
• γglacial: glacial TMI grid

Output

• Aglacial: glacial water-mass matrix

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 mixedlayermask(A,γ)

 function mixedlayermatrix(A, γ, τ)

Read and assemble the circulation matrix from the efficient storage of A and F₀ variables.

Arguments

• A: TMI water-mass matrix
• γ: TMI grid
• τ: uniform residence timescale (years) for all mixed layer points

Output

• Lmix: circulation matrix in xyz format for mixed layer points

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.

oneeastboundary(γ,name=:none,longname="unknown",units="unknown") = ones(1,eastindex(γ),γ,name,longname,units)::BoundaryCondition

onenorthboundary(γ,name=:none,longname="unknown",units="unknown") = ones(2,northindex(γ),γ,name,longname,units)::BoundaryCondition

onesouthboundary(γ,name=:none,longname="unknown",units="unknown") = ones(2,southindex(γ),γ,name,longname,units)::BoundaryCondition

onesurfaceboundary(γ,name=:none,longname="unknown",units="unknown") = ones(3,surfaceindex(γ),γ,name,longname,units)::BoundaryCondition

onewestboundary(γ,name=:none,longname="unknown",units="unknown") = ones(1,westindex(γ),γ,name,longname,units)::BoundaryCondition

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 mat2ncfield

Rename MATLAB variables to NetCDF variables

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 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

update water-mass matrix to agree with new boundaries in γ

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

zeroeastboundary(γ,name=:none,longname="unknown",units="unknown") = zeros(1,eastindex(γ),γ,name,longname,units)::BoundaryCondition

zeronorthboundary(γ,name=:none,longname="unknown",units="unknown") = zeros(2,northindex(γ),γ,name,longname,units)::BoundaryCondition

zerosouthboundary(γ,name=:none,longname="unknown",units="unknown") = zeros(2,southindex(γ),γ,name,longname,units)::BoundaryCondition

zerosurfaceboundary(γ::Grid,name=:none,longname="unknown",units="unknown") = zeros(3,surfaceindex(γ),γ,name,longname,units)::BoundaryCondition

zerowestboundary(γ,name=:none,longname="unknown",units="unknown") = zeros(1,westindex(γ),γ,name,longname,units)::BoundaryCondition

function Γsfc! 
Γsfc! anonymously calls surfacecontrol2field!