DataAssim.jl

x,Hx = FreeRun(ℳ,xi,Q,H,nmax,no)

Performs a free-run with the model ℳ and nmax time-steps starting at the initial condition xi. Observations at the time steps given in no are extracted with the observation operator H.

x,P = KalmanFilter(xi,Pi,ℳ,Q,yo,R,H,nmax,no)

Kalman Filter with the model ℳ and nmax time-steps starting at the initial condition xi and error covariance Pi. Observations yo (and error covariance R) at the time steps given in no are assimilated with the observation operator H.

x,J = fourDVar(
        xi,Pi,ℳ,yo,R,H,nmax,no;
        innerloops = 10,
        outerloops = 2,
        tol = 1e-5)

Incremental 4D-Var with the model ℳ and nmax time-steps starting at the initial condition xi and error covariance Pi with the specified numbers of inner and outer loops. Observations yo (and error covariance R) at the time steps given in no are assimilated with the observation operator H.

Xa,xa = ETKF(Xf,HXf,y,R,H,...)

Computes analysis ensemble Xa based on forecast ensemble Xf and observations y using the ETKF ensemble scheme.

Input arguments:

• Xf: forecast ensemble (n x N)
• HXf: the observation operator applied on the ensemble (product H*Xf)
• y: observations (m)
• R: observation error covariance (m x m).
• H: operator (m x n). Except for the serialEnSRF it is never used and can be empty

Optional keywords arguments:

• debug: set to true to enable debugging. Default (false) is no debugging.
• tolerance: expected rounding error (default 1e-10) for debugging checks. This is not used if debug is false.

Output arguments:

• Xa: the analysis ensemble (n x N)
• xa: the analysis ensemble mean (n)

Notations follows: Sangoma D3.1 http://data-assimilation.net/Documents/sangomaDL3.1.pdf

Xa,xa = EnKF(Xf,HXf,y,R,H,...)

Computes analysis ensemble Xa based on forecast ensemble Xf and observations y using the EnKF ensemble scheme.

Input arguments:

• Xf: forecast ensemble (n x N)
• HXf: the observation operator applied on the ensemble (product H*Xf)
• y: observations (m)
• R: observation error covariance (m x m).
• H: operator (m x n). Except for the serialEnSRF it is never used and can be empty

Optional keywords arguments:

• debug: set to true to enable debugging. Default (false) is no debugging.
• tolerance: expected rounding error (default 1e-10) for debugging checks. This is not used if debug is false.

Output arguments:

• Xa: the analysis ensemble (n x N)
• xa: the analysis ensemble mean (n)

Notations follows: Sangoma D3.1 http://data-assimilation.net/Documents/sangomaDL3.1.pdf

Xa,xa = EnSRF(Xf,HXf,y,R,H,...)

Computes analysis ensemble Xa based on forecast ensemble Xf and observations y using the EnSRF ensemble scheme.

Input arguments:

• Xf: forecast ensemble (n x N)
• HXf: the observation operator applied on the ensemble (product H*Xf)
• y: observations (m)
• R: observation error covariance (m x m).
• H: operator (m x n). Except for the serialEnSRF it is never used and can be empty

Optional keywords arguments:

• debug: set to true to enable debugging. Default (false) is no debugging.
• tolerance: expected rounding error (default 1e-10) for debugging checks. This is not used if debug is false.

Output arguments:

• Xa: the analysis ensemble (n x N)
• xa: the analysis ensemble mean (n)

Notations follows: Sangoma D3.1 http://data-assimilation.net/Documents/sangomaDL3.1.pdf

Xa,xa = EAKF(Xf,HXf,y,R,H,...)

Computes analysis ensemble Xa based on forecast ensemble Xf and observations y using the EAKF ensemble scheme.

Input arguments:

• Xf: forecast ensemble (n x N)
• HXf: the observation operator applied on the ensemble (product H*Xf)
• y: observations (m)
• R: observation error covariance (m x m).
• H: operator (m x n). Except for the serialEnSRF it is never used and can be empty

Optional keywords arguments:

• debug: set to true to enable debugging. Default (false) is no debugging.
• tolerance: expected rounding error (default 1e-10) for debugging checks. This is not used if debug is false.

Output arguments:

• Xa: the analysis ensemble (n x N)
• xa: the analysis ensemble mean (n)

Notations follows: Sangoma D3.1 http://data-assimilation.net/Documents/sangomaDL3.1.pdf

Xa,xa = SEIK(Xf,HXf,y,R,H,...)

Computes analysis ensemble Xa based on forecast ensemble Xf and observations y using the SEIK ensemble scheme.

Input arguments:

• Xf: forecast ensemble (n x N)
• HXf: the observation operator applied on the ensemble (product H*Xf)
• y: observations (m)
• R: observation error covariance (m x m).
• H: operator (m x n). Except for the serialEnSRF it is never used and can be empty

Optional keywords arguments:

• debug: set to true to enable debugging. Default (false) is no debugging.
• tolerance: expected rounding error (default 1e-10) for debugging checks. This is not used if debug is false.

Output arguments:

• Xa: the analysis ensemble (n x N)
• xa: the analysis ensemble mean (n)

Notations follows: Sangoma D3.1 http://data-assimilation.net/Documents/sangomaDL3.1.pdf

Xa,xa = ESTKF(Xf,HXf,y,R,H,...)

Computes analysis ensemble Xa based on forecast ensemble Xf and observations y using the ESTKF ensemble scheme.

Input arguments:

• Xf: forecast ensemble (n x N)
• HXf: the observation operator applied on the ensemble (product H*Xf)
• y: observations (m)
• R: observation error covariance (m x m).
• H: operator (m x n). Except for the serialEnSRF it is never used and can be empty

Optional keywords arguments:

• debug: set to true to enable debugging. Default (false) is no debugging.
• tolerance: expected rounding error (default 1e-10) for debugging checks. This is not used if debug is false.

Output arguments:

• Xa: the analysis ensemble (n x N)
• xa: the analysis ensemble mean (n)

Notations follows: Sangoma D3.1 http://data-assimilation.net/Documents/sangomaDL3.1.pdf

Xa,xa = serialEnSRF(Xf,HXf,y,R,H,...)

Computes analysis ensemble Xa based on forecast ensemble Xf and observations y using the serialEnSRF ensemble scheme.

Input arguments:

• Xf: forecast ensemble (n x N)
• HXf: the observation operator applied on the ensemble (product H*Xf)
• y: observations (m)
• R: observation error covariance (m x m).
• H: operator (m x n). Except for the serialEnSRF it is never used and can be empty

Optional keywords arguments:

• debug: set to true to enable debugging. Default (false) is no debugging.
• tolerance: expected rounding error (default 1e-10) for debugging checks. This is not used if debug is false.

Output arguments:

• Xa: the analysis ensemble (n x N)
• xa: the analysis ensemble mean (n)

Notations follows: Sangoma D3.1 http://data-assimilation.net/Documents/sangomaDL3.1.pdf

Xa,xa = local_ETKF(Xf,H,y,diagR,part,selectObs,...)

Computes analysis ensemble Xa based on forecast ensemble Xf using the observation y using the local ETKF.

Inputs:

• Xf: forecast ensemble (n x N)
• H: observation operator (m x n)
• y: observation (m x 1)
• diagR: diagonal of the observation error covariance R (m x 1)
• part: vector of integer "labels". Every element of the state vector with the same number belong to the same subdomain
• selectObs: callback routine to select observations with a within a subdomain. As input is takes an integer representing the index of the state vector and returns a vector of weights (m x 1). For example:

     selectObs(i) = exp(- ((x[i] - xobs[:]).^2 + (y(i) - yobs[:]).^2)/L^2 );

or

     selectObs(i) = compact_locfun(L,...
         sqrt((x[i] - xobs[:]).^2 + (y[i] - yobs[:]).^2));

where x and y is the horizontal model grid, xobs and yobs are the locations of the observations and L is a correlation length-scale

Optional inputs:

• display: if true, then display progress (false is the default)
• minweight: analysis is performed using observations for which weights is larger than minweight. (default 1e-8)
• HXf: if non empty, then it is the product H*Xf. In this case, H is not used

Output:

• Xa: the analysis ensemble (n x N)
• xa`: the analysis ensemble mean (n x 1)

See also: compact_locfun

Xa,xa = local_EnKF(Xf,H,y,diagR,part,selectObs,...)

Computes analysis ensemble Xa based on forecast ensemble Xf using the observation y using the local EnKF.

Inputs:

• Xf: forecast ensemble (n x N)
• H: observation operator (m x n)
• y: observation (m x 1)
• diagR: diagonal of the observation error covariance R (m x 1)
• part: vector of integer "labels". Every element of the state vector with the same number belong to the same subdomain
• selectObs: callback routine to select observations with a within a subdomain. As input is takes an integer representing the index of the state vector and returns a vector of weights (m x 1). For example:

     selectObs(i) = exp(- ((x[i] - xobs[:]).^2 + (y(i) - yobs[:]).^2)/L^2 );

or

     selectObs(i) = compact_locfun(L,...
         sqrt((x[i] - xobs[:]).^2 + (y[i] - yobs[:]).^2));

where x and y is the horizontal model grid, xobs and yobs are the locations of the observations and L is a correlation length-scale

Optional inputs:

• display: if true, then display progress (false is the default)
• minweight: analysis is performed using observations for which weights is larger than minweight. (default 1e-8)
• HXf: if non empty, then it is the product H*Xf. In this case, H is not used

Output:

• Xa: the analysis ensemble (n x N)
• xa`: the analysis ensemble mean (n x 1)

See also: compact_locfun

Xa,xa = local_EnSRF(Xf,H,y,diagR,part,selectObs,...)

Computes analysis ensemble Xa based on forecast ensemble Xf using the observation y using the local EnSRF.

Inputs:

• Xf: forecast ensemble (n x N)
• H: observation operator (m x n)
• y: observation (m x 1)
• diagR: diagonal of the observation error covariance R (m x 1)
• part: vector of integer "labels". Every element of the state vector with the same number belong to the same subdomain
• selectObs: callback routine to select observations with a within a subdomain. As input is takes an integer representing the index of the state vector and returns a vector of weights (m x 1). For example:

     selectObs(i) = exp(- ((x[i] - xobs[:]).^2 + (y(i) - yobs[:]).^2)/L^2 );

or

     selectObs(i) = compact_locfun(L,...
         sqrt((x[i] - xobs[:]).^2 + (y[i] - yobs[:]).^2));

where x and y is the horizontal model grid, xobs and yobs are the locations of the observations and L is a correlation length-scale

Optional inputs:

• display: if true, then display progress (false is the default)
• minweight: analysis is performed using observations for which weights is larger than minweight. (default 1e-8)
• HXf: if non empty, then it is the product H*Xf. In this case, H is not used

Output:

• Xa: the analysis ensemble (n x N)
• xa`: the analysis ensemble mean (n x 1)

See also: compact_locfun

Xa,xa = local_EAKF(Xf,H,y,diagR,part,selectObs,...)

Computes analysis ensemble Xa based on forecast ensemble Xf using the observation y using the local EAKF.

Inputs:

• Xf: forecast ensemble (n x N)
• H: observation operator (m x n)
• y: observation (m x 1)
• diagR: diagonal of the observation error covariance R (m x 1)
• part: vector of integer "labels". Every element of the state vector with the same number belong to the same subdomain
• selectObs: callback routine to select observations with a within a subdomain. As input is takes an integer representing the index of the state vector and returns a vector of weights (m x 1). For example:

     selectObs(i) = exp(- ((x[i] - xobs[:]).^2 + (y(i) - yobs[:]).^2)/L^2 );

or

     selectObs(i) = compact_locfun(L,...
         sqrt((x[i] - xobs[:]).^2 + (y[i] - yobs[:]).^2));

where x and y is the horizontal model grid, xobs and yobs are the locations of the observations and L is a correlation length-scale

Optional inputs:

• display: if true, then display progress (false is the default)
• minweight: analysis is performed using observations for which weights is larger than minweight. (default 1e-8)
• HXf: if non empty, then it is the product H*Xf. In this case, H is not used

Output:

• Xa: the analysis ensemble (n x N)
• xa`: the analysis ensemble mean (n x 1)

See also: compact_locfun

Xa,xa = local_SEIK(Xf,H,y,diagR,part,selectObs,...)

Computes analysis ensemble Xa based on forecast ensemble Xf using the observation y using the local SEIK.

Inputs:

• Xf: forecast ensemble (n x N)
• H: observation operator (m x n)
• y: observation (m x 1)
• diagR: diagonal of the observation error covariance R (m x 1)
• part: vector of integer "labels". Every element of the state vector with the same number belong to the same subdomain
• selectObs: callback routine to select observations with a within a subdomain. As input is takes an integer representing the index of the state vector and returns a vector of weights (m x 1). For example:

     selectObs(i) = exp(- ((x[i] - xobs[:]).^2 + (y(i) - yobs[:]).^2)/L^2 );

or

     selectObs(i) = compact_locfun(L,...
         sqrt((x[i] - xobs[:]).^2 + (y[i] - yobs[:]).^2));

where x and y is the horizontal model grid, xobs and yobs are the locations of the observations and L is a correlation length-scale

Optional inputs:

• display: if true, then display progress (false is the default)
• minweight: analysis is performed using observations for which weights is larger than minweight. (default 1e-8)
• HXf: if non empty, then it is the product H*Xf. In this case, H is not used

Output:

• Xa: the analysis ensemble (n x N)
• xa`: the analysis ensemble mean (n x 1)

See also: compact_locfun

Xa,xa = local_ESTKF(Xf,H,y,diagR,part,selectObs,...)

Computes analysis ensemble Xa based on forecast ensemble Xf using the observation y using the local ESTKF.

Inputs:

• Xf: forecast ensemble (n x N)
• H: observation operator (m x n)
• y: observation (m x 1)
• diagR: diagonal of the observation error covariance R (m x 1)
• part: vector of integer "labels". Every element of the state vector with the same number belong to the same subdomain
• selectObs: callback routine to select observations with a within a subdomain. As input is takes an integer representing the index of the state vector and returns a vector of weights (m x 1). For example:

     selectObs(i) = exp(- ((x[i] - xobs[:]).^2 + (y(i) - yobs[:]).^2)/L^2 );

or

     selectObs(i) = compact_locfun(L,...
         sqrt((x[i] - xobs[:]).^2 + (y[i] - yobs[:]).^2));

where x and y is the horizontal model grid, xobs and yobs are the locations of the observations and L is a correlation length-scale

Optional inputs:

• display: if true, then display progress (false is the default)
• minweight: analysis is performed using observations for which weights is larger than minweight. (default 1e-8)
• HXf: if non empty, then it is the product H*Xf. In this case, H is not used

Output:

• Xa: the analysis ensemble (n x N)
• xa`: the analysis ensemble mean (n x 1)

See also: compact_locfun

Abstract base-class of models. A model should implement forecast step, tangent-linear and adjoint step

ℳ = ModelMatrix(M)

Linear model defined by the matrix M.

ℳ = ModelFun(nonlinear_forecast,tangent_linear_model,adjoint_model)

Model defined by the functions nonlinear_forecast,tangent_linear_model and adjoint_model.

ℳ = LinShallowWater1DModel(dt,g,h,L,imax)

Linear 1D shallow water model.

Example

dt = 1.
g = 9.81
h = 100
imax = 101
L = 10000
LinShallowWater1DModel(dt,g,h,L,imax)

ℳ = Lorenz63Model(dt,σ=10.,β = 8/3.,ρ = 28.)

Lorenz, 1963 model[1] integrated with a 2nd order Runge-Kutta scheme.

[1] https://doi.org/10.1175/1520-0469(1963)020%3C0130:DNF%3E2.0.CO;2

 fun = compact_locfun(r)

Smooth compact localization function at the (scaled) distance r. fun is zero if r > 2 and one if r is 0. (Gaspari et al. (1999), equation 4.10, [1])

[1] http://dx.doi.org/10.1002/qj.49712555417