Utilities

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 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 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 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 ones(dim::Int64,dimval::Int64,γ::Grid)::BoundaryCondition

   Initialize boundary condition with ones

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 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 *(c::Field,d::Field)::Field`
Field by field multiplication is element-by-element.

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.

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

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

function location_obs(field, locs, γ)

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