osse/EnKF/EnKF.F

c***  ensemble kalman filter subroutine

        subroutine EnKF(A,z,H,ns,no,varmes,nrsamp)

c***  A is the matrix of ensemble forecasts, z is the observation,
c***  H is the observation operator, ns is the size of the model
c***  state, no is the number of observation locations, varmes is
c***  the observational uncertainty, and nrsamp is the ensemble size.

        integer   ns, no, nrsamp
        integer   ijim

        real, intent(inout)         :: A(ns,nrsamp)
        real, intent(in)            :: z(no)

        integer i,j,m,mm,info,idum,ibawk
        real    rcond
        real    H(no,ns)
        real    vave(ns), vvar(ns), varmes(no)
        real    temp

        real, allocatable    :: K(:,:)
        real, allocatable    :: RR(:,:)
        real, allocatable    :: RRR(:,:)
        real, allocatable    :: w(:)
        real, allocatable    :: tmp(:,:)
        real, allocatable    :: d(:)
        real, allocatable    :: pd(:,:)
        real, allocatable    :: transA(:,:)
        real, allocatable    :: transtmp(:,:)
        real, allocatable    :: tmp2(:,:)


        allocate (K(no,ns), RR(no,no), RRR(no,no))
        allocate (transA(nrsamp,ns),transtmp(nrsamp,no))
        allocate (w(ns), tmp(no,nrsamp))
        allocate (d(no), pd(1,ns), tmp2(1,ns))

c***  initialising the observational error covariance matrix and
c***  the projector from model space to observational space.
         RR=0.0

c***  define obs error covariance matrix as a diagonal matrix
         do m=1,no
          RR(m,m)=varmes(m)
         enddo

c***  calculate ensemble mean
         call ranmean2(A,w,ns,nrsamp)

c***  remove ensemble mean from the ensemble members
         do i=1,nrsamp
          A(:,i)=A(:,i)-w(:)
         enddo

c***  tmp takes on the product of the observation projector and 
c***  the ensemble anomalies.
         tmp=matmul(H,A)

c***  K takes on the value of H(A-ABAR)(A-ABAR)^T
         K=matmul(tmp,transpose(A))

c***  K is what is called B in eqn (6)
         K=(1.0/(float(nrsamp)-1.0))*K

c***  restore A by adding back ensemble mean value.
         do i=1,nrsamp
          A(:,i)=A(:,i)+w(:)
         enddo  

c***  RR is obs error cov matrix, and the forecast error cov is added to
c***  it, making it the first half of the right hand side of equation (8).
c***  if one was to wipe out spurious long distance correlations, it would
c***  happen here.
        RR=RR+(1./(float(nrsamp)-1.))*matmul(tmp,transpose(tmp))

c***  standard solver for coefficients beta
         call spoco(RR,no,no,rcond,w,info)

c***  correct each member of the ensemble
         do j=1,nrsamp

          call random2(d,no)
          d=z+sqrt(varmes)*d-matmul(H,A(:,j))
          call sposl(RR,no,no,d)
          A(:,j)=A(:,j)+matmul(d,K) 

         enddo                  !A has now been corrected by the data.

        deallocate (K, RR, RRR)
        deallocate (transA, transtmp)
        deallocate (w, tmp)
        deallocate (d,pd,tmp2)

        return
        end


1	c*** ensemble kalman filter subroutine
2
3	subroutine EnKF(A,z,H,ns,no,varmes,nrsamp)
4
5	c*** A is the matrix of ensemble forecasts, z is the observation,
6	c*** H is the observation operator, ns is the size of the model
7	c*** state, no is the number of observation locations, varmes is
8	c*** the observational uncertainty, and nrsamp is the ensemble size.
9
10	integer ns, no, nrsamp
11	integer ijim
12
13	real, intent(inout) :: A(ns,nrsamp)
14	real, intent(in) :: z(no)
15
16	integer i,j,m,mm,info,idum,ibawk
17	real rcond
18	real H(no,ns)
19	real vave(ns), vvar(ns), varmes(no)
20	real temp
21
22	real, allocatable :: K(:,:)
23	real, allocatable :: RR(:,:)
24	real, allocatable :: RRR(:,:)
25	real, allocatable :: w(:)
26	real, allocatable :: tmp(:,:)
27	real, allocatable :: d(:)
28	real, allocatable :: pd(:,:)
29	real, allocatable :: transA(:,:)
30	real, allocatable :: transtmp(:,:)
31	real, allocatable :: tmp2(:,:)
32
33
34	allocate (K(no,ns), RR(no,no), RRR(no,no))
35	allocate (transA(nrsamp,ns),transtmp(nrsamp,no))
36	allocate (w(ns), tmp(no,nrsamp))
37	allocate (d(no), pd(1,ns), tmp2(1,ns))
38
39	c*** initialising the observational error covariance matrix and
40	c*** the projector from model space to observational space.
41	RR=0.0
42
43	c*** define obs error covariance matrix as a diagonal matrix
44	do m=1,no
45	RR(m,m)=varmes(m)
46	enddo
47
48	c*** calculate ensemble mean
49	call ranmean2(A,w,ns,nrsamp)
50
51	c*** remove ensemble mean from the ensemble members
52	do i=1,nrsamp
53	A(:,i)=A(:,i)-w(:)
54	enddo
55
56	c*** tmp takes on the product of the observation projector and
57	c*** the ensemble anomalies.
58	tmp=matmul(H,A)
59
60	c*** K takes on the value of H(A-ABAR)(A-ABAR)^T
61	K=matmul(tmp,transpose(A))
62
63	c*** K is what is called B in eqn (6)
64	K=(1.0/(float(nrsamp)-1.0))*K
65
66	c*** restore A by adding back ensemble mean value.
67	do i=1,nrsamp
68	A(:,i)=A(:,i)+w(:)
69	enddo
70
71	c*** RR is obs error cov matrix, and the forecast error cov is added to
72	c*** it, making it the first half of the right hand side of equation (8).
73	c*** if one was to wipe out spurious long distance correlations, it would
74	c*** happen here.
75	RR=RR+(1./(float(nrsamp)-1.))*matmul(tmp,transpose(tmp))
76
77	c*** standard solver for coefficients beta
78	call spoco(RR,no,no,rcond,w,info)
79
80	c*** correct each member of the ensemble
81	do j=1,nrsamp
82
83	call random2(d,no)
84	d=z+sqrt(varmes)*d-matmul(H,A(:,j))
85	call sposl(RR,no,no,d)
86	A(:,j)=A(:,j)+matmul(d,K)
87
88	enddo !A has now been corrected by the data.
89
90	deallocate (K, RR, RRR)
91	deallocate (transA, transtmp)
92	deallocate (w, tmp)
93	deallocate (d,pd,tmp2)
94
95	return
96	end
97
98
99
100
101