| 1 |
afe |
1.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 |
|
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph