Output: the correction vector du0
n = len(u0b)
#determine the assimilation window
t = t[:ind_m[-1]+1] #cut the time till the last observation point
nt = len(t)-1
dt = t[1] - t[0]
ub = np.zeros([n,nt+1]) #base trajectory
Ri = np.linalg.inv(R)
ub[:,0] = u0b
U = np.eye(n,n) #Initialization of U
Q = np.zeros((1,n)) #Dh*U
ef = np.zeros((1,1)) #w-h(u)
W = np.zeros((1,1)) #weighting matrix
km = 0
nt_m = len(ind_m)
if opt == 0: #Euler
#forward model
for k in range(nt):
ub[:,k+1] = euler(rhs,ub[:,k],dt,*args)
DM = Jeuler(rhs,Jrhs,ub[:,k],dt,*args)
U = DM @ U
if (km<nt_m) and (k+1==ind_m[km]):
tmp = w[:,km] - ObsOp(ub[:,k+1])
ek = tmp.reshape(-1,1)
```
*Fluids* **2020**, *5*, 225

```
ef = np.vstack((ef,ek))
Qk = JObsOp(ub[:,k+1]) @ U
Q = np.vstack((Q,Qk))
W = block_diag(W,Ri)
km = km + 1
elif opt == 1: #RK4