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 write(unit=stdout,fmt=*) ' '
51
52 !--------------------------------------------------------------------------
53 ! [2.0]: Test global eigenvectors:
54 !--------------------------------------------------------------------------
55
56 ! [2.1] Print out global eigenvectors:
57
58 write(unit=stdout,fmt='(2A)') '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 write(unit=stdout,fmt='(A)') " "
65
66 ! [2.2]: Test eigenvector orthogonality: sum_k (e_m(k) e_n(k)) = delta_mn:
67
68 allocate(array_mask(1:kz))
69 allocate(matrix2(1:kz,1:kz))
70
71 write(unit=stdout,fmt='(2A)') &
72 'Eigenvector orthogonality check for ', trim(name)
73 write(unit=stdout,fmt='(A)')' Mode Diagonal Maximum off-diagonal'
74 do k1 = 1, kz
75 do k2 = 1, kz
76 matrix2(k1,k2) = sum(be_eigenvec(:,k1) * be_eigenvec(:,k2))
77 end do
78
79 array_mask(1:kz) =.true.
80 array_mask(k1) = .false.
81 max_off_diag = MAXVAL(ABS(matrix2(k1,:)),mask=array_mask(:))
82 write(unit=stdout,fmt='(I4,4x,1pe12.4,10x,1pe12.4)')k1, matrix2(k1,k1), &
83 max_off_diag
84 end do
85 write(unit=stdout,fmt=*) ' '
86
87 ! [2.3] Test eigenvectors completeness - sum_m (e_m(k1) e_m(k2)) = delta_k1k2
88
89 write(unit=stdout,fmt='(2A)') &
90 'Eigenvector completeness check for ', trim(name)
91 write(unit=stdout,fmt='(A)')' Level Diagonal Maximum off-diagonal'
92 do k1 = 1, kz
93 do k2 = 1, kz
94 matrix2(k1,k2) = sum(be_eigenvec(k1,:) * be_eigenvec(k2,:))
95 end do
96
97 array_mask(1:kz) =.true.
98 array_mask(k1) = .false.
99 max_off_diag = MAXVAL(ABS(matrix2(k1,:)),mask=array_mask(:))
100 write(unit=stdout,fmt='(I4,4x,1pe12.4,10x,1pe12.4)')k1, matrix2(k1,k1), &
101 max_off_diag
102 end do
103 write(unit=stdout,fmt=*) ' '
104
105 !-------------------------------------------------------------------------
106 ! [3.0]: Tidy up:
107 !-------------------------------------------------------------------------
108
109 deallocate(matrix2)
110 deallocate(array_mask)
111
112 if (trace_use) call da_trace_exit("da_check_eof_decomposition")
113
114 end subroutine da_check_eof_decomposition
115
116