da_check_eof_decomposition.inc
References to this file elsewhere.
1 subroutine da_check_eof_decomposition(be_eigenval, be_eigenvec, name)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Check eigenvectors E of vertical covariance matrix B_{x} which
5 ! have been read in from NMC-statistics file.
6 !
7 ! Method: E and L (eigenvalues) are computed using LAPACK/BLAS software in
8 ! the NMC code using the definition E^{T} B_{x} E = L. This routine checks
9 ! for eigenvector orthogonality and completeness as defined below.
10 !---------------------------------------------------------------------------
11
12 implicit none
13
14 real, intent(in) :: be_eigenval(:) ! Back. error eigenvalues.
15 real, intent(in) :: be_eigenvec(:,:) ! Back. error eigenvector
16 character(len=*), intent(in) :: name ! Variable name.
17
18 integer :: kz ! Size of 3rd dimension.
19 integer :: k, k1, k2, m ! Loop counters
20 real :: tot_variance ! Total variance.
21 real :: cum_variance ! Cumulative variance.
22 real :: max_off_diag ! Maximum off-diagonal.
23
24 real, allocatable :: matrix2(:,:) ! 2D Work matrix.
25 logical, allocatable :: array_mask(:) ! Array mask for MAXVAL.
26
27 if (trace_use) call da_trace_entry("da_check_eof_decomposition")
28
29 !----------------------------------------------------------------------
30 ! [1.0]: Initialise:
31 !----------------------------------------------------------------------
32
33 kz = size(be_eigenval)
34
35 !----------------------------------------------------------------------
36 ! [2.0]: Print out global eigenvalues (used for truncation):
37 !----------------------------------------------------------------------
38
39 cum_variance = 0.0
40 tot_variance = sum(be_eigenval(1:kz))
41
42 write(unit=stdout,fmt='(A)')' Name Mode Eigenvalue Normalised Variance'
43 do k = 1, kz
44 cum_variance = be_eigenval(k) + cum_variance
45 write(unit=stdout,fmt='(A,I4,e14.4,f10.4)') &
46 trim(name), k, be_eigenval(k), cum_variance / tot_variance
47 end do
48 write(unit=stdout,fmt=*) ' '
49 write(unit=stdout,fmt='(A,e13.5)')' Total variance = Tr(Bv) = ', tot_variance
50
51 !--------------------------------------------------------------------------
52 ! [2.0]: Test global eigenvectors:
53 !--------------------------------------------------------------------------
54
55 ! [2.1] Print out global eigenvectors:
56
57 write(unit=stdout,fmt='(3A)') ' da_check_eof_decomposition:', &
58 ' Domain eigenvectors for ', trim(name)
59
60 write(unit=stdout,fmt='(50i13)')(m, m=1,kz)
61 do k = 1, kz
62 write(unit=stdout,fmt='(I3,50e13.5)')k, (be_eigenvec(k,m), m=1,kz)
63 end do
64
65 ! [2.2]: Test eigenvector orthogonality: sum_k (e_m(k) e_n(k)) = delta_mn:
66
67 allocate(array_mask(1:kz))
68 allocate(matrix2(1:kz,1:kz))
69
70 write(unit=stdout,fmt='(2A)') &
71 ' Eigenvector orthogonality check for ', trim(name)
72 write(unit=stdout,fmt='(A)')' Mode Diagonal Maximum off-diagonal'
73 do k1 = 1, kz
74 do k2 = 1, kz
75 matrix2(k1,k2) = sum(be_eigenvec(:,k1) * be_eigenvec(:,k2))
76 end do
77
78 array_mask(1:kz) =.true.
79 array_mask(k1) = .false.
80 max_off_diag = MAXVAL(ABS(matrix2(k1,:)),mask=array_mask(:))
81 write(unit=stdout,fmt='(I4,4x,1pe12.4,10x,1pe12.4)')k1, matrix2(k1,k1), &
82 max_off_diag
83 end do
84 write(unit=stdout,fmt=*) ' '
85
86 ! [2.3] Test eigenvectors completeness - sum_m (e_m(k1) e_m(k2)) = delta_k1k2
87
88 write(unit=stdout,fmt='(2A)') &
89 ' Eigenvector completeness check for ', trim(name)
90 write(unit=stdout,fmt='(A)')' Level Diagonal Maximum off-diagonal'
91 do k1 = 1, kz
92 do k2 = 1, kz
93 matrix2(k1,k2) = sum(be_eigenvec(k1,:) * be_eigenvec(k2,:))
94 end do
95
96 array_mask(1:kz) =.true.
97 array_mask(k1) = .false.
98 max_off_diag = MAXVAL(ABS(matrix2(k1,:)),mask=array_mask(:))
99 write(unit=stdout,fmt='(I4,4x,1pe12.4,10x,1pe12.4)')k1, matrix2(k1,k1), &
100 max_off_diag
101 end do
102 write(unit=stdout,fmt=*) ' '
103
104 !-------------------------------------------------------------------------
105 ! [3.0]: Tidy up:
106 !-------------------------------------------------------------------------
107
108 deallocate(matrix2)
109 deallocate(array_mask)
110
111 if (trace_use) call da_trace_exit("da_check_eof_decomposition")
112
113 end subroutine da_check_eof_decomposition
114
115