Compute a variational analysis of arbitrarily located observations.
[fi,err,s] = divand(mask,pmn,xi,x,f,len,lambda,...);
[fi,err] = divand(mask,pmn,xi,x,f,len,lambda,...);
[fi,s] = divand(mask,pmn,xi,x,f,len,lambda,...);
Perform an n-dimensional variational analysis of the observations f located at
the coordinates x. The array fi represent the interpolated field at the grid
defined by the coordinates xi and the scales factors pmn.
Input:
mask: binary mask delimiting the domain. 1 is inside and 0 outside.
For oceanographic application, this is the land-sea mask.
pmn: scale factor of the grid. pmn is a cell array with n elements. Every
element represents the scale factor of the corresponding dimension. Its
inverse is the local resolution of the grid in a particular dimension.
xi: cell array with n elements. Every element represents a coordinate
of the final grid on which the observations are interpolated
x: cell array with n elements. Every element represents a coordinate of
the observations
f: value of the observations *minus* the background estimate (m-by-1 array).
(see note)
len: correlation length
lambda: signal-to-noise ratio of observations (if lambda is a scalar).
The larger this value is, the closer is the field fi to the
observation.
if lambda is a scalar:
R = 1/lambda I, where R is the observation error covariance
matrix),
if lambda is a vector
a vector (R = diag(lambda)) or a matrix or a CovarParam object (R =
lambda).
Optional input arguments specified as pairs of keyword and values:
'velocity', vel: velocity of advection constraint. The default is
no-advection constraint
'alpha': alpha is vector of coefficients multiplying various terms in the
cost function. The first element multiplies the norm.
The other i-th element of alpha multiplies the (i+1)-th derivative.
Per default, the highest derivative is m = ceil(1+n/2) where n is the
dimension of the problem.
The values of alpha is the (m+1)th row of the Pascal triangle:
m=0 1
m=1 1 1
m=1 1 2 1 (n=1,2)
m=2 1 3 3 1 (n=3,4)
...
'diagnostics': 0 or 1 turns diagnostic and debugging information on (1) or
off (0, default). If on, they will be returned as the last output
argument
'EOF', EOF: sub-space constraint. Orthogonal (EOF' WE^2 EOF = I) (units of
EOF: m^(-n/2))
'EOF_scaling', EOF_scaling: (dimensional)
'constraint': a structure with user specified constrain
'moddim': modulo for cyclic dimension (vector with n elements).
Zero is used for non-cyclic dimensions. Halo points should
not be included for cyclic dimensions. For example if the first dimension
is cyclic, then the grid point corresponding to mask(1,j) should be
between mask(end,1) (left neighbor) and mask(2,j) (right neighbor)
'fracdim': fractional indices (n-by-m array). If this array is specified,
then x and xi are not used.
'inversion': direct solver ('chol' for Cholesky factorization) or a
interative solver ('pcg' for preconditioned conjugate gradient) can be
used.
'compPC': function that returns a preconditioner for the primal formulation
if inversion is set to 'pcg'. The function has the following arguments:
[M1,M2] = compPC(iB,H,R)
where iB is the inverse background error covariance, H the observation
operator and R the error covariance of the observation. The used
preconditioner M will be M = M1 * M2 (or just M = M1 if M2 is empty).
Per default a modified incomplete Cholesky factorization will be used a
preconditioner.
Note: 'velocity' and 'constraint' may appear multiple times
Output:
fi: the analysed field
err: error variance of the analysis field relative to the error variance of
the background
s: structure
s.iB: adimensional
s.E: scaled EOF (dimensional)
Note:
If zero is not a valid first guess for your variable (as it is the case for
e.g. ocean temperature), you have to subtract the first guess from the
observations before calling divand and then add the first guess back in.
Example:
see divand_simple_example.m
Package: divand
