da_etkf.inc
References to this file elsewhere.
1 subroutine da_etkf(ndim,nanals,nobs,ens,ens_ob,oberrvar,obs,nout,naccumt1, &
2 naccumt2,nstartaccum1,nstartaccum2,tainflatinput,rhoinput)
3
4 !-----------------------------------------------------------------------
5 ! Purpose: subroutine for ETKF perturbation updates
6 ! Xuguang Wang, Jan. 2006
7 !
8 ! references:
9 ! Bishop et al. 2001 MWR
10 ! Wang and Bishop 2003 JAS
11 ! Wang et al. 2004 MWR
12 ! Wang et al. 2005 MWR
13 ! 1) nanals, ensemble size
14 ! 2) ndim, dimension of the perturbation vector that needs to be updated
15 ! 3) nobs, number of observations assmiilated
16 ! 4) ens, array of perturbations before (Xf') and after (Xa') ETKF update
17 ! 5) ens_ob, array of HXf
18 ! 6) oberrvar, observation error variance, listed in the same sequence as HXf
19 ! 7) obs, observations assmiilated
20 ! 8) naccumt1, number of previous cycles immediately before the current cycle,
21 ! which is needed for calculating adaptive inflation factor. naccumt1 < 0 for
22 ! pre-specified inflation
23 ! 9) naccumt2, number of previous cycles immediately before the current cycle,
24 ! which is needed for calculating the rho factor in the latest version of
25 ! ETKF. naccumt2 < 0 for using the older version of the ETKF. naccumt2 = 0
26 ! for using pre-specified rho factor
27 ! 10) nstartaccum1, the cycle from which the accumulation of previous
28 ! naccumt1 cycles in 8) starts
29 ! 11) nstartaccum2, the cycle from which the accumulation of previous naccumt2
30 ! in 9) starts
31 ! 12) tainflatinput, pre-specified inflation, if not using adaptive inflation
32 ! 13) rhoinput, pre-specified rho factor, if not using adaptively determined
33 ! rho factor
34 ! 14) nout, record number for output squared innovations and the forecast
35 ! error variance projected onto ensemble subspace, which is related to 8) and 9)
36 !-----------------------------------------------------------------------
37
38 implicit none
39
40 integer, intent(in) :: nanals,ndim,nobs
41 real, intent(inout), dimension(ndim,nanals) :: ens
42 real, intent(inout), dimension(nobs,nanals) :: ens_ob
43 real, intent(in), dimension(nobs) :: oberrvar
44 real, intent(in), dimension(nobs) :: obs
45 real, intent(in) :: tainflatinput,rhoinput
46 integer, intent(in) :: nout,naccumt1,naccumt2,nstartaccum1,nstartaccum2
47
48 real, dimension(ndim) :: ensmean
49 real, dimension(nobs) :: ensmean_ob
50 real, dimension(nobs) :: obsinc
51 integer :: ij, nanal, i, j, k, nt
52 integer :: info, lwork
53 real, allocatable, dimension(:) :: work
54 real, dimension(1) :: workt
55 real, dimension(nanals) :: eignv, eignv1
56 real, dimension(nanals, nanals) :: hzthz, C, TR
57 real, dimension(nanals-1, nanals) :: CT
58 real, dimension(nanals, nanals-1) :: cgamma, T
59 real, dimension(nobs,nanals-1) :: E
60 real, dimension(nobs,nanals) :: enspert_ob_tmp
61 ! file for output eigen values filename_eignv
62 ! file for output factors filename_factor
63 character (len=150) :: filename_factor,filename_eignv
64 ! file for squared innovation vectors
65 ! file for squared projection of fcst err variance onto ensemble subspace
66 character (len=150) :: filename_inno2,filename_proj2
67 real, dimension(nanals-1) :: proj
68 real :: ainflat, tainflat, rho
69 real :: proj2, proj2sum, proj2sum1
70 real :: squareinnosum, squareinno, squareinnosum1
71 real :: ta, tb, tracehpfht, aftersum_eignv, sum_eignv
72 filename_factor = "factor.dat"
73 filename_eignv = "eignv.dat"
74 filename_inno2 = "inno2.dat"
75 filename_proj2 = "proj2.dat"
76
77 ! open file for output factors ainflat, tainflat, rho
78 open(109,form="unformatted",file=filename_factor,access="direct",recl=4)
79
80 ! open file for output eigenvalues
81 open(119,form="unformatted",file=filename_eignv,access="direct",recl=4*nanals)
82
83 ! open file for
84 open(99,form="unformatted",file=filename_inno2,access="direct",recl=4)
85
86 ! open file for
87 open(89,form="unformatted",file=filename_proj2,access="direct",recl=4)
88
89 call cpu_time(ta)
90
91 ! compute ensemble mean.
92 do ij=1,ndim
93 ensmean(ij) = sum(ens(ij,:))/float(nanals)
94 end do
95
96 do ij=1,nobs
97 ensmean_ob(ij) = sum(ens_ob(ij,:))/float(nanals)
98 end do
99
100 ! remove ens mean to get ensemble perturbations
101
102 do nanal=1,nanals
103 ens(:,nanal) = ens(:,nanal)-ensmean(:)
104 ens_ob(:,nanal) = ens_ob(:,nanal)-ensmean_ob(:)
105 end do
106
107 ! calculate HZTHZ Wang and Bishop 2003
108 do i = 1, nobs
109 enspert_ob_tmp(i,:) = ens_ob(i,:)/sqrt(oberrvar(i))
110 end do
111
112 call innerprod(enspert_ob_tmp,hzthz,nobs,nanals)
113 hzthz = hzthz/float(nanals-1)
114
115 ! calculate C and Gamma in Wang and Bishop 2003
116 ! in outputs, hzthz contains C, eignv contains gamma
117 call SSYEV('V', 'L', nanals , hzthz, nanals, eignv1, workt, -1, info)
118 lwork = int(workt(1))
119 allocate (work(lwork))
120 call SSYEV('V', 'L', nanals , hzthz, nanals, eignv1, work, lwork, info)
121 deallocate(work)
122
123 ! note eignv1 output from SSYEV is in ascending order !!!
124 ! re-order both eigenvalues and eigenvectors first
125
126 do i = 1, nanals
127 eignv(i) = eignv1(nanals-i+1)
128 end do
129
130 write(119,rec=nout) eignv
131
132 write (unit=stdout,fmt=*) "eignv=",eignv
133
134 ! note eigenvectors output from SSYEV are stored in columns !!!
135 do i = 1, nanals
136 C(:,i) = hzthz(:,nanals-i+1)
137 end do
138
139
140 ! calculate the inflation factor used to inflate the ETKF initial perts
141 ! by ainflat = (obsinc*obsinc/oberrvar-nobs)/trace(pfht_ob/oberrvar),
142 ! tainflate = ainflate1*ainflat2*.... in Wang and Bishop 2003
143
144 tracehpfht = 0.0
145 do i = 1, nanals-1
146 tracehpfht = tracehpfht + eignv(i)
147 end do
148
149 obsinc = ensmean_ob - obs
150
151 squareinno = 0.0
152 do i = 1, nobs
153 squareinno = squareinno + obsinc(i)*obsinc(i)/oberrvar(i)
154 end do
155
156 write(99,rec=nout) squareinno
157
158 write (unit=stdout,fmt=*) "squareinno=", squareinno
159 write (unit=stdout,fmt=*) "nout=", nout
160 write (unit=stdout,fmt=*) "nstartaccum1=", nstartaccum1
161 write (unit=stdout,fmt=*) "tracehpfht=",tracehpfht
162 write (unit=stdout,fmt=*) "nobs=",nobs
163
164 if (naccumt1 .gt. 0) then
165 squareinnosum = squareinno
166
167 if (nout.ge.nstartaccum1) then
168 write (unit=stdout,fmt=*) "OK"
169 squareinnosum = 0.0
170 do nt = nout-naccumt1+1, nout
171 read(99,rec=nt) squareinnosum1
172 write (unit=stdout,fmt=*) "squareinnosum1=", squareinnosum1
173 squareinnosum = squareinnosum + squareinnosum1/float(naccumt1)
174 end do
175 end if
176
177 write (unit=stdout,fmt=*) "squareinnosum=",squareinnosum
178
179 ainflat = (squareinnosum-float(nobs))/tracehpfht
180 if (ainflat .le. 0) then
181 ainflat = 1.0
182 end if
183
184 if (nout .le. nstartaccum1) then
185 tainflat = ainflat
186 end if
187
188 if (nout .gt. nstartaccum1) then
189 read(109, rec=(nout-2)*3+2) tainflat
190 tainflat = tainflat*ainflat
191 end if
192
193 write(109, rec=(nout-1)*3+1) ainflat
194 write(109, rec=(nout-1)*3+2) tainflat
195
196 write (unit=stdout,fmt=*) "tainflat=", tainflat
197 write (unit=stdout,fmt=*) "ainflat=", ainflat
198
199 else
200 ! this is for pre-specified inflation factor
201 tainflat = tainflatinput
202 write (unit=stdout,fmt=*) "tainflat=", tainflat
203 write(109, rec=(nout-1)*3+2) tainflat
204 end if
205
206
207 ! calculate the rho factor following Wang et al 2005
208 if (naccumt2 .gt. 0) then
209 ! calculate E, the ensemble subspace eign vectors in Wang and Bishop (2003)
210 ! E=R(-1/2)HZCgamma^(-1/2)
211
212 do i = 1, nanals
213 do j = 1, nanals-1
214 cgamma(i,j) = C(i,j)*sqrt(1.0/eignv(j))
215 end do
216 end do
217
218 ! R(-1/2)HZ
219 do i = 1, nobs
220 enspert_ob_tmp(i,:) = ens_ob(i,:)/sqrt(oberrvar(i))/sqrt(float(nanals-1))
221 end do
222
223 call matmulti(enspert_ob_tmp,cgamma,E,nobs,nanals-1,nanals)
224
225 ! project normalized (divided by oberrstdev) innovation vector onto E
226 proj = 0.0
227 do i = 1, nanals-1
228 do k = 1, nobs
229 proj(i) = proj(i) + obsinc(k)/sqrt(oberrvar(k))*E(k,i)
230 end do
231 end do
232
233 ! get rho = (sum(proj*proj)-dim. of proj))/(normalized innovation^2-nobs)
234 ! since nanals is relatively small, need follow wang et al (2005) for accumulation
235 ! relative error = sqrt(2/(dim of proj)) = e.g., sqrt(2/207)=10%
236 ! e.g., 50mem 2wks = sqrt(2/(14*49)) = 5% = e.g. 3wks sqrt(2/(21*49))=4%
237
238 proj2 = sum(proj*proj)
239
240 write(89,rec=nout) proj2
241
242 proj2sum = proj2
243
244 if (nout .ge. nstartaccum2) then
245 proj2sum = 0.0
246 do nt = nout-naccumt2+1,nout
247 read(89,rec=nt) proj2sum1
248 proj2sum = proj2sum + proj2sum1/float(naccumt2)
249 end do
250 end if
251
252 squareinnosum = squareinno
253
254 if (nout.ge.nstartaccum2) then
255 write (unit=stdout,fmt=*) "OK"
256 squareinnosum = 0.0
257 do nt = nout-naccumt2+1, nout
258 read(99,rec=nt) squareinnosum1
259 write (unit=stdout,fmt=*) "squareinnosum1=", squareinnosum1
260 squareinnosum = squareinnosum + squareinnosum1/float(naccumt2)
261 end do
262 end if
263
264 write (unit=stdout,fmt=*) "squareinnosum=",squareinnosum
265
266 rho = (proj2sum-(float(nanals-1)))/(squareinnosum-float(nobs))
267 if (rho .le. 0) then
268 rho = 1.0
269 end if
270
271 write(109, rec=(nout-1)*3+3) rho
272
273 else if (naccumt2 .lt. 0.0) then
274 ! this is for the old formulation of ETKF in Wang and Bishop 2003
275 rho = 1.0
276 else
277 ! this if for pre-specified rho factor
278 rho = rhoinput
279 end if
280
281 write (unit=stdout,fmt=*) "rho=",rho
282 write (unit=stdout,fmt=*) "proj2sum=",proj2sum
283
284 ! calculate the grand transformation matrix following Wang et al 2005
285 ! TR = C*(rho*gamma+I)^(-1/2)*C^T
286 do i = 1, nanals
287 do j = 1, nanals-1
288 T(i,j) = C(i,j)*sqrt(1.0/(rho*eignv(j)+1.0))
289 end do
290 end do
291
292 do i = 1, nanals-1
293 do j = 1, nanals
294 CT(i,j) = C(j,i)
295 end do
296 end do
297
298 call matmulti(T,CT,TR,nanals,nanals,nanals-1)
299 ! apply ainflat
300 TR = sqrt(tainflat)*TR
301
302 write (unit=stdout,fmt=*) "ens"
303 write (unit=stdout,fmt=*) maxval(ens)
304 write (unit=stdout,fmt=*) minval(ens)
305 write (unit=stdout,fmt=*) "eigen/(eigen*rho+1)"
306 write (unit=stdout,fmt=*) eignv/(eignv*rho+1.0)
307 sum_eignv = 0.0
308 aftersum_eignv = 0.0
309 do i = 1, nanals
310 sum_eignv = sum_eignv + eignv(i)
311 aftersum_eignv = aftersum_eignv + eignv(i)/(eignv(i)*rho+1.0)
312 end do
313 write (unit=stdout,fmt=*) "sum_eignv=", sum_eignv
314 write (unit=stdout,fmt=*) "sum_aftereignv=", aftersum_eignv
315
316
317 write (unit=stdout,fmt=*) "avg_ens=",sum(ens)/float(ndim*nanals)
318 write (unit=stdout,fmt=*) "std_ens=",sum(ens*ens)/float(ndim*nanals)
319
320 ! calculate the ETKF initial perturbations
321 call matmultiover(ens,TR,ndim,nanals)
322 write (unit=stdout,fmt=*) "ens"
323 write (unit=stdout,fmt=*) maxval(ens)
324 write (unit=stdout,fmt=*) minval(ens)
325 write (unit=stdout,fmt=*) "afteravg_ens=",sum(ens)/float(ndim*nanals)
326 write (unit=stdout,fmt=*) "afterstd_ens=",sum(ens*ens)/float(ndim*nanals)
327
328 ! calculate the ETKF initial perturbations
329
330 call cpu_time(tb)
331 write (unit=stdout,fmt=*) "total cpu time =", tb-ta
332 write (unit=stdout,fmt=*) "total cpu time tb =", tb
333
334 end subroutine da_etkf
335
336