Main Structure of the Filter Algorithm

3.2.5 Resampling

Since the EnKF updates in the analysis phase the whole ensemble of model states the algorithm can proceed directly to the next ensemble forecast without the need of a re-sampling algorithm. In contrast to this, a new state ensemble representing ˇP^a_k and x^a_k has to be generated when the SEIK filter is used. This can be done by a transforma-tion of the forecast ensemble. Applying the SEEK filter, the forecasted modes of the covariance matrix can be used directly in the next forecast phase. In general, these are no more the basis vectors of the error subspace, since they are not orthonormal. A re-orthonormalization of the modes is recommendable and can be performed occasion-ally to stabilize the mode forecast. The choice whether an algorithm with or without re-initialization is used has no particular implications for the performance of the filter algorithms.

Subroutine SEIK Main(n,N,X) int n {state dimension, input} int N {ensemble size, input}

real X(n, N) {state ensemble array, input} real x(n) {state estimate}

real Uinv(N −1, N −1) {inverse of eigenvalue matrix} int i {ensemble loop counter}

int step {time step counter}

int m {dimension of observation vector} real t_a {physical time}

1: call User Analysis seik(0,n,N,X) {call to user analysis routine}

2: loop

3: call Next Observation(step, nsteps, t_a)

{get number of time steps, user supplied}

4: if nsteps= 0 then

5: exit loop

6: end if

7: for i=1 to N do

8: call Interface Evolver(n,X(N),nsteps, t_a) {forecast state vector, user supplied}

9: end for

10: step←step+nsteps

11: call User Analysis(−step,n,N,X) {call to user supplied analysis routine}

12: call SEIK Analysis(step,n,N,x,Uinv,X) {perform filter analysis phase}

13: call SEIK Resample(n,N,x,Uinv,X) {perform ensemble resampling}

14: call User Analysis(step,n,N,X) {call to user supplied analysis routine}

15: end loop

Algorithm 3.1: Structure of the filter main subroutine for the SEIK algorithm. The arraysxand Uinvare required for the resampling computed inSEIK Resample. They are initialized in the analysis routine SEIK Analysis.

The calls to the subroutine User Analysis in algorithm 3.1 provide the possibility to examine the assimilation progress during the execution. Here the user can analyze either the forecast or the analysis state ensemble. To distinguish both cases, the sub-routine is called with the negative of the time step indexstepsin the forecast case. The routine permits, e.g., to compute ensemble means or variances estimated by the filter.

In addition, the ensemble or analysis quantities can be written to files. For physical consistency it can be necessary to post-process the analysis states, for example to en-sure mass conservation of a model ocean. This post-processing can be also performed in User Analysis when called after the filter analysis phase.

In the forecast phase an ensemble ofN model state vectorsX ={⁽¹⁾x, . . . ,^(N)x}is evolved for nsteps time steps from the model time t_a to the time t_b =t_a+nsteps·∆t where ∆tis the time step size. This requires to perform N model evolutions beginning

from the same model time t_a. The ensemble forecast is controlled by the filter, since the model does not need to consider filter details. The parameters nsteps and t_a are dependent on the data assimilation problem rather than on the model or the filter algorithm. Thus, they have to be provided by the user. For flexibility and to achieve a clear structure we implement the initialization of nsteps and t_a by a call to the user supplied subroutine Next Observation. It has as input the current time step step.

Outputs are nsteps and t_a.

Having obtained the values of nsteps and t_a, the forecast is performed in a loop over all ensemble vectors. Each of the vectors is handed over to the subroutine Inter-face Evolver together with the stepping information. This interface routine initializes the state fields of the model from the state vector and calls the time stepper routine of the model. Finally the fields are written back into the state vector and the routine returns. Since Interface Evolver is model dependent it has to be supplied by the user.

The forecast phase requires that the N model evolutions are independent. Thus, any reused variables of the model have to be re-initialized.

Subsequent to the forecast phase, the analysis will be computed. In algorithm 3.1 this is performed in the subroutine SEIK Analysis. We discuss the implementation of the analysis phases of the three filters in the following section. Finally, the ensemble will be resampled in the SEIK algorithm. The new ensemble is computed in the subroutine SEIK Resample. The implementation of the resampling phases of SEIK and SEEK is described in section 3.3.3.

Subroutine SEEK Main(n,r,x,Uinv,V) int r {rank, input}

real x(n) {state estimate, input}

real Uinv(r, r) {inverse of eigenvalue matrix, input} real V(n, r) {mode matrix, input}

real ² {coefficient for gradient approximation} ...

1: for i=1 tor do

2: V(:, r)←x+²V(:, r) {generate ensemble from modes}

3: end for

4: for i=1 tor do

5: call Interface Evolver(n,V(:,r),nsteps, t_a) {forecast ensemble vector}

6: end for

7: call Interface Evolver(n,x,nsteps, t_a) {forecast central state vector}

8: for i=1 tor do

9: V(:, r)←²⁻¹(x−V(:, r) {generate forecast modes from ensemble}

10: end for ...

Algorithm 3.2: Structure of forecast part of the filter main subroutine for the SEEK algorithm

The structure of the main routine of the EnKF algorithm is analogous to that of the SEIK filter and thus not shown. The only functional difference is that the EnKF algorithm does not call a resampling routine. Further, the arraysUinv and x are not required. For the SEEK algorithm the forecast part is different from the two other algorithms. In SEEK the state vector x and the mode matrix V are evolved. The structure of the forecast loop using a gradient approximation for the evolution of the modes stored in V is shown in algorithm 3.2.