gen_be_stage2_1dvar.f90
References to this file elsewhere.
1 program gen_be_stage2_1dvar
2
3 !---------------------------------------------------------------------------------------
4 ! Purpose: Calculate multivariate covariances required in 1D-Var via "NMC-method" or
5 ! ensemble approaches.
6 !
7 ! Method: Accumulates statistics in loop over times/ensemble members.
8 ! History:
9 ! Creation (from gen_be_stage_2) Dale Barker 01/07/2005.
10 !---------------------------------------------------------------------------------------
11
12 use da_control, only : filename_len
13 use da_gen_be, only : da_eof_decomposition, da_eof_decomposition_test
14 use da_tools_serial, only : da_get_unit,da_advance_cymdh
15
16 implicit none
17
18 real, parameter :: t_factor = 1.0 ! Normalization factor.
19 real, parameter :: q_factor = 1000.0 ! Normalization factor.
20 real, parameter :: ps_factor = 0.01 ! Normalization factor.
21 real, parameter :: u_factor = 1.0 ! Normalization factor.
22
23 character*10 :: start_date, end_date ! Starting and ending dates.
24 character*10 :: date, new_date ! Current date (ccyymmddhh).
25 character*10 :: variable ! Variable name
26 character(len=filename_len) :: dat_dir ! Input data directory.
27 character*80 :: expt ! Experiment ID.
28 character(len=filename_len) :: filename ! Input filename.
29 character*3 :: ce ! Member index -> character.
30 integer :: ni, nj, nk, nkdum ! Grid dimensions.
31 integer :: nk1 ! nk + 1.
32 integer :: nkbe ! 2*nk1 + 3.
33 integer :: i, j, k, member, k2 ! Loop counters.
34 integer :: kbe, k2be ! Loop counters.
35 integer :: b ! Bin marker.
36 integer :: sdate, cdate, edate ! Starting, current ending dates.
37 integer :: interval ! Period between dates (hours).
38 integer :: ne ! Number of ensemble members.
39 integer :: bin_type ! Type of bin to average over.
40 integer :: num_bins ! Number of bins (3D fields).
41 integer :: num_bins2d ! Number of bins (2D fields).
42 real :: lat_min, lat_max ! Used if bin_type = 2 (degrees).
43 real :: binwidth_lat ! Used if bin_type = 2 (degrees).
44 real :: hgt_min, hgt_max ! Used if bin_type = 2 (m).
45 real :: binwidth_hgt ! Used if bin_type = 2 (m).
46 real :: coeffa, coeffb ! Accumulating mean coefficients.
47 logical :: first_time ! True if first file.
48 logical :: ldum1, ldum2 ! Dummy logicals.
49 logical :: testing_eofs ! Test EOF calculation in inverse if true.
50
51 real, allocatable :: temp(:,:,:) ! Temperature.
52 real, allocatable :: q(:,:,:) ! Humidity.
53 real, allocatable :: ps(:,:) ! Surface pressure.
54 real, allocatable :: u10(:,:) ! 10m wind.
55 real, allocatable :: v10(:,:) ! 10m wind.
56 integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point.
57 integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point.
58 integer, allocatable:: bin_pts2d(:) ! Number of points in bin (2D fields).
59 real, allocatable :: covar_tt(:,:,:) ! Covariance between input fields.
60 real, allocatable :: covar_tq(:,:,:) ! Covariance between input fields.
61 real, allocatable :: covar_tps(:,:) ! Covariance between input fields.
62 real, allocatable :: covar_tu10(:,:) ! Covariance between input fields.
63 real, allocatable :: covar_tv10(:,:) ! Covariance between input fields.
64 real, allocatable :: covar_qq(:,:,:) ! Covariance between input fields.
65 real, allocatable :: covar_qps(:,:) ! Covariance between input fields.
66 real, allocatable :: covar_qu10(:,:) ! Covariance between input fields.
67 real, allocatable :: covar_qv10(:,:) ! Covariance between input fields.
68 real, allocatable :: covar_psps(:) ! Covariance between input fields.
69 real, allocatable :: covar_psu10(:) ! Covariance between input fields.
70 real, allocatable :: covar_psv10(:) ! Covariance between input fields.
71 real, allocatable :: covar_u10u10(:) ! Covariance between input fields.
72 real, allocatable :: covar_u10v10(:) ! Covariance between input fields.
73 real, allocatable :: covar_v10v10(:) ! Covariance between input fields.
74 real, allocatable :: be(:,:,:) ! BE matrix.
75 real, allocatable :: be_inv(:,:,:) ! BE matrix inverse.
76 real, allocatable :: be_loc(:,:) ! BE matrix.
77 real, allocatable :: be_loc_inv(:,:) ! BE matrix inverse.
78 real, allocatable :: be_test(:,:) ! BE matrix inverse.
79
80 namelist / gen_be_stage2_1dvar_nl / start_date, end_date, interval, &
81 ne, testing_eofs, expt, dat_dir
82
83 integer :: ounit,iunit,namelist_unit
84
85 !---------------------------------------------------------------------------------------------
86 write(6,'(a)')' [1] Initialize namelist variables and other scalars.'
87 !---------------------------------------------------------------------------------------------
88
89 call da_get_unit(ounit)
90 call da_get_unit(iunit)
91 call da_get_unit(namelist_unit)
92
93 start_date = '2004030312'
94 end_date = '2004033112'
95 interval = 24
96 ne = 1
97 testing_eofs = .false.
98 expt = 'gen_be_stage2_1dvar'
99 dat_dir = '/mmmtmp1/dmbarker'
100
101 open(unit=namelist_unit, file='gen_be_stage2_1dvar_nl.nl', &
102 form='formatted', status='old', action='read')
103 read(namelist_unit, gen_be_stage2_1dvar_nl)
104 close(namelist_unit)
105
106 read(start_date(1:10), fmt='(i10)')sdate
107 read(end_date(1:10), fmt='(i10)')edate
108 write(6,'(4a)')' Computing multivariate covariances for 1D-Var'
109 write(6,'(4a)') ' Time period is ', start_date, ' to ', end_date
110 write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
111 write(6,'(a,i8)')' Number of ensemble members at each time = ', ne
112
113 date = start_date
114 cdate = sdate
115
116 !---------------------------------------------------------------------------------------------
117 write(6,'(2a)')' [2] Read fields, and calculate correlations'
118 !---------------------------------------------------------------------------------------------
119
120 write(6,'(/a)')' Unbiased perturbations read in have been scaled by the following factors:'
121 write(6,'(a,f15.5)')' T scaling factor = ', t_factor
122 write(6,'(a,f15.5)')' q scaling factor = ', q_factor
123 write(6,'(a,f15.5)')' ps scaling factor = ', ps_factor
124 write(6,'(a,f15.5)')' u scaling factor = ', u_factor
125
126 first_time = .true.
127
128 do while ( cdate <= edate )
129 write(6,'(a,a)')' Processing data for date ', date
130
131 do member = 1, ne
132
133 write(ce,'(i3.3)')member
134
135 ! Read T:
136 variable = 't'
137 filename = trim(variable)//'/'//date(1:10)
138 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
139 open (iunit, file = filename, form='unformatted')
140 read(iunit)ni, nj, nk
141 nk1 = nk + 1
142
143 if ( first_time ) then
144 write(6,'(a,3i8)')' i, j, k dimensions are ', ni, nj, nk1
145 write(6,'(a,3i8)')' Note k dimension expanded by one to include T2m, q2m in array.'
146
147 allocate( bin(1:ni,1:nj,1:nk) ) ! Just needed for read in, not used.
148 allocate( temp(1:ni,1:nj,1:nk1) )
149 allocate( q(1:ni,1:nj,1:nk1) )
150 allocate( bin2d(1:ni,1:nj) )
151 allocate( ps(1:ni,1:nj) )
152 allocate( u10(1:ni,1:nj) )
153 allocate( v10(1:ni,1:nj) )
154 end if
155
156 read(iunit)temp(1:ni,1:nj,2:nk+1)
157 close(iunit)
158
159 if ( first_time ) then
160 ! Read bin info:
161 filename = 'bin.data'
162 open (iunit, file = filename, form='unformatted')
163 read(iunit)bin_type
164 read(iunit)lat_min, lat_max, binwidth_lat
165 read(iunit)hgt_min, hgt_max, binwidth_hgt
166 read(iunit)num_bins, num_bins2d
167 read(iunit)bin(1:ni,1:nj,1:nk) ! Just read in, not used.
168 read(iunit)bin2d(1:ni,1:nj)
169 close(iunit)
170
171 allocate( bin_pts2d(1:num_bins2d) )
172 allocate( covar_tt(1:nk1,1:nk1,1:num_bins2d) ) ! Only store one quadrant of matrix.
173 allocate( covar_tq(1:nk1,1:nk1,1:num_bins2d) ) ! The only nk1 * nk1 matrix
174 allocate( covar_tps(1:nk1,1:num_bins2d) )
175 allocate( covar_tu10(1:nk1,1:num_bins2d) )
176 allocate( covar_tv10(1:nk1,1:num_bins2d) )
177 allocate( covar_qq(1:nk1,1:nk1,1:num_bins2d) ) ! Only store one quadrant of matrix.
178 allocate( covar_qps(1:nk1,1:num_bins2d) )
179 allocate( covar_qu10(1:nk1,1:num_bins2d) )
180 allocate( covar_qv10(1:nk1,1:num_bins2d) )
181 allocate( covar_psps(1:num_bins2d) )
182 allocate( covar_psu10(1:num_bins2d) )
183 allocate( covar_psv10(1:num_bins2d) )
184 allocate( covar_u10u10(1:num_bins2d) )
185 allocate( covar_u10v10(1:num_bins2d) )
186 allocate( covar_v10v10(1:num_bins2d) )
187 bin_pts2d(:) = 0
188 covar_tt(:,:,:) = 0.0
189 covar_tq(:,:,:) = 0.0
190 covar_tps(:,:) = 0.0
191 covar_tu10(:,:) = 0.0
192 covar_tv10(:,:) = 0.0
193 covar_qq(:,:,:) = 0.0
194 covar_qps(:,:) = 0.0
195 covar_qu10(:,:) = 0.0
196 covar_qv10(:,:) = 0.0
197 covar_psps(:) = 0.0
198 covar_psu10(:) = 0.0
199 covar_psv10(:) = 0.0
200 covar_u10u10(:) = 0.0
201 covar_u10v10(:) = 0.0
202 covar_v10v10(:) = 0.0
203 first_time = .false.
204 end if
205
206 ! Read t2:
207 variable = 't2'
208 filename = trim(variable)//'/'//date(1:10)
209 filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
210 open (iunit, file = filename, form='unformatted')
211 read(iunit)ni, nj, nkdum
212 read(iunit)ldum1, ldum2 ! Dummy logicals.
213 read(iunit)temp(1:ni,1:nj,1) ! Read 2m temperature into k=1 of expanded temp array.
214 close(iunit)
215
216 ! Read q:
217 variable = 'q'
218 filename = trim(variable)//'/'//date(1:10)
219 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
220 open (iunit, file = filename, form='unformatted')
221 read(iunit)ni, nj, nk
222 read(iunit)q(1:ni,1:nj,2:nk+1)
223 close(iunit)
224
225 ! Read q2:
226 variable = 'q2'
227 filename = trim(variable)//'/'//date(1:10)
228 filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
229 open (iunit, file = filename, form='unformatted')
230 read(iunit)ni, nj, nkdum
231 read(iunit)ldum1, ldum2 ! Dummy logicals.
232 read(iunit)q(1:ni,1:nj,1) ! Read 2m humidity into k=1 of expanded q array.
233 close(iunit)
234
235 ! Read ps:
236 variable = 'ps'
237 filename = trim(variable)//'/'//date(1:10)
238 filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
239 open (iunit, file = filename, form='unformatted')
240 read(iunit)ni, nj, nkdum
241 read(iunit)ldum1, ldum2 ! Dummy logicals.
242 read(iunit)ps
243 close(iunit)
244
245 ! Read u10:
246 variable = 'u10'
247 filename = trim(variable)//'/'//date(1:10)
248 filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
249 open (iunit, file = filename, form='unformatted')
250 read(iunit)ni, nj, nkdum
251 read(iunit)ldum1, ldum2 ! Dummy logicals.
252 read(iunit)u10
253 close(iunit)
254
255 ! Read v10:
256 variable = 'v10'
257 filename = trim(variable)//'/'//date(1:10)
258 filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
259 open (iunit, file = filename, form='unformatted')
260 read(iunit)ni, nj, nkdum
261 read(iunit)ldum1, ldum2 ! Dummy logicals.
262 read(iunit)v10
263 close(iunit)
264
265 ! Calculate multivariate covariances:
266
267 do j = 1, nj
268 do i = 1, ni
269 b = bin2d(i,j)
270 bin_pts2d(b) = bin_pts2d(b) + 1
271 coeffa = 1.0 / real(bin_pts2d(b))
272 coeffb = real(bin_pts2d(b)-1) * coeffa
273
274 do k = 1, nk1
275
276 ! T/T covariance (symmetric):
277 do k2 = 1, k
278 covar_tt(k,k2,b) = coeffb * covar_tt(k,k2,b) + &
279 coeffa * temp(i,j,k) * temp(i,j,k2)
280 end do
281
282 ! T/q covariance:
283 do k2 = 1, nk1
284 covar_tq(k,k2,b) = coeffb * covar_tq(k,k2,b) + &
285 coeffa * temp(i,j,k) * q(i,j,k2)
286 end do
287
288 ! T/ps covariance:
289 covar_tps(k,b) = coeffb * covar_tps(k,b) + coeffa * temp(i,j,k) * ps(i,j)
290
291 ! T/u10 covariance:
292 covar_tu10(k,b) = coeffb * covar_tu10(k,b) + coeffa * temp(i,j,k) * u10(i,j)
293
294 ! T/u10 covariance:
295 covar_tv10(k,b) = coeffb * covar_tv10(k,b) + coeffa * temp(i,j,k) * v10(i,j)
296
297 ! q/q covariance (symmetric):
298 do k2 = 1, k
299 covar_qq(k,k2,b) = coeffb * covar_qq(k,k2,b) + &
300 coeffa * q(i,j,k) * q(i,j,k2)
301 end do
302
303 ! q/ps covariance:
304 covar_qps(k,b) = coeffb * covar_qps(k,b) + coeffa * q(i,j,k) * ps(i,j)
305
306 ! q/u10 covariance:
307 covar_qu10(k,b) = coeffb * covar_qu10(k,b) + coeffa * q(i,j,k) * u10(i,j)
308
309 ! q/v10 covariance:
310 covar_qv10(k,b) = coeffb * covar_qv10(k,b) + coeffa * q(i,j,k) * v10(i,j)
311 end do
312
313 ! 1D variances:
314 covar_psps(b) = coeffb * covar_psps(b) + coeffa * ps(i,j) * ps(i,j)
315 covar_psu10(b) = coeffb * covar_psu10(b) + coeffa * ps(i,j) * u10(i,j)
316 covar_psv10(b) = coeffb * covar_psv10(b) + coeffa * ps(i,j) * v10(i,j)
317 covar_u10u10(b) = coeffb * covar_u10u10(b) + coeffa * u10(i,j) * u10(i,j)
318 covar_u10v10(b) = coeffb * covar_u10v10(b) + coeffa * u10(i,j) * v10(i,j)
319 covar_v10v10(b) = coeffb * covar_v10v10(b) + coeffa * v10(i,j) * v10(i,j)
320
321 end do
322 end do
323 end do ! End loop over ensemble members.
324
325 ! Calculate next date:
326 call da_advance_cymdh( date, interval, new_date )
327 date = new_date
328 read(date(1:10), fmt='(i10)')cdate
329 end do ! End loop over times.
330
331 !---------------------------------------------------------------------------------------------
332 write(6,'(2a)')' [3] Pack into total background error covariance and write '
333 !---------------------------------------------------------------------------------------------
334
335 nkbe = 2 * nk1 + 3
336 allocate( be(1:nkbe,1:nkbe,1:num_bins2d) )
337 be = 0.0
338
339 do b = 1, num_bins2d
340
341 ! Pack TT:
342 do k = 1, nk1
343 kbe = k
344 do k2 = 1, k
345 k2be = k2
346 be(kbe,k2be,b) = covar_tt(k,k2,b)
347 end do
348 end do
349
350 ! Pack Tq:
351 do k = 1, nk1
352 kbe = nk1 + k
353 do k2 = 1, nk1
354 k2be = k2
355 be(kbe,k2be,b) = covar_tq(k,k2,b)
356 end do
357 end do
358
359 ! Pack Tps:
360 kbe = 2 * nk1 + 1
361 do k2 = 1, nk1
362 k2be = k2
363 be(kbe,k2be,b) = covar_tps(k2,b)
364 end do
365
366 ! Pack Tu10:
367 kbe = 2 * nk1 + 2
368 do k2 = 1, nk1
369 k2be = k2
370 be(kbe,k2be,b) = covar_tu10(k2,b)
371 end do
372
373 ! Pack Tv10:
374 kbe = 2 * nk1 + 3
375 do k2 = 1, nk1
376 k2be = k2
377 be(kbe,k2be,b) = covar_tv10(k2,b)
378 end do
379
380 ! Pack qq:
381 do k = 1, nk1
382 kbe = nk1 + k
383 do k2 = 1, k
384 k2be = nk1 + k2
385 be(kbe,k2be,b) = covar_qq(k,k2,b)
386 end do
387 end do
388
389 ! Pack qps:
390 kbe = 2 * nk1 + 1
391 do k2 = 1, nk1
392 k2be = nk1 + k2
393 be(kbe,k2be,b) = covar_qps(k2,b)
394 end do
395
396 ! Pack qu10:
397 kbe = 2 * nk1 + 2
398 do k2 = 1, nk1
399 k2be = nk1 + k2
400 be(kbe,k2be,b) = covar_qu10(k2,b)
401 end do
402
403 ! Pack qv10:
404 kbe = 2 * nk1 + 3
405 do k2 = 1, nk1
406 k2be = nk1 + k2
407 be(kbe,k2be,b) = covar_qv10(k2,b)
408 end do
409
410 ! Pack psps:
411 kbe = 2 * nk1 + 1
412 k2be = 2 * nk1 + 1
413 be(kbe,k2be,b) = covar_psps(b)
414
415 ! Pack psu10:
416 kbe = 2 * nk1 + 2
417 k2be = 2 * nk1 + 1
418 be(kbe,k2be,b) = covar_psu10(b)
419
420 ! Pack psv10:
421 kbe = 2 * nk1 + 3
422 k2be = 2 * nk1 + 1
423 be(kbe,k2be,b) = covar_psv10(b)
424
425 ! Pack u10u10:
426 kbe = 2 * nk1 + 2
427 k2be = 2 * nk1 + 2
428 be(kbe,k2be,b) = covar_u10u10(b)
429
430 ! Pack u10v10:
431 kbe = 2 * nk1 + 3
432 k2be = 2 * nk1 + 2
433 be(kbe,k2be,b) = covar_u10v10(b)
434
435 ! Pack v10v10:
436 kbe = 2 * nk1 + 3
437 k2be = 2 * nk1 + 3
438 be(kbe,k2be,b) = covar_v10v10(b)
439
440 ! Fill in auto-covariances by symmetry:
441 do k = 1, nkbe
442 do k2 = k+1, nkbe
443 be(k,k2,b) = be(k2,k,b)
444 end do
445 end do
446
447 end do
448
449 ! Write data for input in 1D-Var:
450 filename = 'gen_be_stage2_1dvar.dat'
451 open (ounit, file = filename, form='unformatted')
452 write(ounit)ni, nj, nk1
453 write(ounit)num_bins2d
454 write(ounit)be
455 close(ounit)
456
457 !---------------------------------------------------------------------------------------------
458 write(6,'(2a)')' [4] Calculate and write inverse covariances '
459 !---------------------------------------------------------------------------------------------
460
461 allocate( be_inv(1:nkbe,1:nkbe,1:num_bins2d) )
462 allocate( be_loc(1:nkbe,1:nkbe) )
463 allocate( be_loc_inv(1:nkbe,1:nkbe) )
464 allocate( be_test(1:nkbe,1:nkbe) )
465 be_test = 0.0
466
467 do b = 1, num_bins2d
468
469 be_loc(:,:) = be(:,:,b)
470 call da_get_be_inverse( testing_eofs, nkbe, be_loc, be_loc_inv )
471 be_inv(:,:,b) = be_loc_inv(:,:)
472
473 ! Check inverse:
474 be_test = matmul(be_loc,be_loc_inv)
475
476 ! Output test covariance matrix (should be identity matrix):
477 do k = 1, nkbe
478 do k2 = 1, nkbe
479 write(52,'(f15.5)')be_test(k,k2)
480 end do
481 end do
482 end do
483
484 filename = 'gen_be_stage2_1dvar.inv.dat'
485 open (ounit, file = filename, form='unformatted')
486 write(ounit)ni, nj, nk1
487 write(ounit)num_bins2d
488 write(ounit)be_inv
489 close(ounit)
490
491 !---------------------------------------------------------------------------------------------
492 write(6,'(2a)')' [4] Read in again, and output for graphics '
493 !---------------------------------------------------------------------------------------------
494
495 write(6,'(a,i10)')' num_bins2d = ', num_bins2d
496
497 be = 0.0
498 be_inv = 0.0
499
500 ! Output BE covariance matrix:
501 filename = 'gen_be_stage2_1dvar.dat'
502 open (iunit, file = filename, form='unformatted', status = 'old')
503 read(iunit)ni, nj, nk1
504 read(iunit)num_bins2d
505 read(iunit)be
506 close(iunit)
507
508 do b = 1, num_bins2d
509 do k = 1, nkbe
510 do k2 = 1, nkbe
511 write(50,'(f15.5)')be(k,k2,b)
512 end do
513 end do
514 end do
515
516 ! Output inverse BE covariance matrix:
517 filename = 'gen_be_stage2_1dvar.inv.dat'
518 open (iunit, file = filename, form='unformatted', status = 'old')
519 read(iunit)ni, nj, nk1
520 read(iunit)num_bins2d
521 read(iunit)be_inv
522 close(iunit)
523
524 do b = 1, num_bins2d
525 do k = 1, nkbe
526 do k2 = 1, nkbe
527 write(51,'(f15.5)')be_inv(k,k2,b)
528 end do
529 end do
530 end do
531
532 contains
533
534 subroutine da_get_be_inverse( testing_eofs, nk, be, be_inv )
535
536 ! Get inverse of a square matrix.
537 ! Dale Barker: 2005/11/03.
538
539 implicit none
540
541 real, parameter :: variance_threshold = 1e-6 ! Percentage of variance discarded.
542
543 logical, intent(inout) :: testing_eofs ! Test if true.
544 integer, intent(in) :: nk ! Array dimension.
545 real, intent(in) :: be(1:nk,1:nk) ! Input square matrix.
546 real, intent(out) :: be_inv(1:nk,1:nk) ! Output inverse matrix.
547
548 integer :: m, k, k2 ! Loop counters.
549 integer :: mmax ! Maximum mode (after variance truncation).
550
551 real :: summ ! Summation dummy.
552 real :: total_var_inv ! 1 / Total variance of matrix.
553 real :: cumul_variance ! Cumulative variance of matrix.
554
555 real :: evec(1:nk,1:nk) ! Eigenvectors.
556 real :: eval(1:nk) ! Eigenvalues.
557 real :: eval_inv(1:nk) ! 1 / Eigenvalues.
558 real :: LamInvET(1:nk,1:nk) ! ET/sqrt(Eigenvalue).
559
560 ! Calculate eigenvectors/values:
561 call da_eof_decomposition( nk, be, evec, eval )
562
563 ! Optionally test calculation:
564 if ( testing_eofs ) then
565 call da_eof_decomposition_test( nk, be, evec, eval )
566 testing_eofs = .false.
567 end if
568
569 ! Truncate eigenvalues to ensure inverse is not dominated by rounding error:
570 summ = 0.0
571 do m = 1, nk
572 summ = summ + eval(m)
573 end do
574 total_var_inv = 1.0 / summ
575
576 cumul_variance = 0.0
577 mmax = nk
578 do m = 1, nk
579 cumul_variance = cumul_variance + eval(m) * total_var_inv
580 if ( eval(m) < 0.0 .or. cumul_variance > 1.0 - variance_threshold ) then
581 mmax = m - 1
582 exit
583 end if
584 end do
585
586 ! L{-1} . E^T:
587 LamInvET(:,:) = 0.0
588 eval_inv(1:mmax) = 1.0 / eval(1:mmax)
589 do k = 1, nk
590 do m = 1, mmax
591 LamInvET(m,k) = evec(k,m) * eval_inv(m)
592 end do
593 end do
594
595 ! B^-1 = E . L{-1} . E^T:
596 do k = 1, nk
597 do k2 = 1, k
598 summ = 0.0
599 do m = 1, nk
600 summ = summ + evec(k,m) * LamInvET(m,k2)
601 end do
602 be_inv(k,k2) = summ
603 end do
604 end do
605
606 ! Fill in rest of matrix by symmetry:
607 do k = 1, nk
608 do k2 = k+1, nk ! Symmetry.
609 be_inv(k,k2) = be_inv(k2,k)
610 end do
611 end do
612
613 end subroutine da_get_be_inverse
614
615 end program gen_be_stage2_1dvar