da_eof_decomposition_test.inc
References to this file elsewhere.
1 subroutine da_eof_decomposition_test (kz, bx, e, l)
2
3 !------------------------------------------------------------------------------
4 ! Purpose:
5 ! [1] Print eigenvalues:
6 ! [2] Test orthogonality of eigenvectors - sum_k (e_m(k) e_n(k)) = delta_mn:
7 ! [3] Test eigenvectors completeness - sum_m (e_m(k1) e_m(k2)) = delta_k1k2:
8 ! [4] Check B correctness: B = E*L*E^T
9 !------------------------------------------------------------------------------
10
11 implicit none
12
13 integer, intent(in) :: kz ! Dimension of BE matrix
14 real, intent(in) :: bx(1:kz,1:kz) ! Global vert. background error.
15 real, intent(in) :: e(1:kz,1:kz) ! Eigenvectors of Bx.
16 real, intent(in) :: l(1:kz) ! Eigenvalues of Bx.
17
18 integer :: k, k1, k2, m ! Loop counters
19 real :: tot_variance ! Total variance.
20 real :: cum_variance ! Cumulative variance.
21 real :: max_off_diag ! Maximum off-diagonal.
22
23 real :: work(1:kz,1:kz) ! 2D Work matrix.
24 real :: bc(1:kz,1:kz) ! 2D Work matrix.
25 logical :: array_mask(1:kz) ! Array mask for MAXVAL.
26
27 if (trace_use) call da_trace_entry("da_eof_decomposition_test")
28
29 !-------------------------------------------------------------------------
30 ! [1] Print eigenvalues:
31 !-------------------------------------------------------------------------
32
33 tot_variance = sum(l(1:kz))
34 cum_variance = 0.0
35
36 write(unit=stdout,fmt='(A)')' Mode Eigenvalue Cumulative Variance e(k,k)'
37
38 do k = 1, kz
39 cum_variance = cum_variance + l(k)
40 write(unit=stdout,fmt='(I4,4x,e12.4,10x,f8.4,4x,e12.4)') &
41 k, l(k), cum_variance / tot_variance, e(k,k)
42 end do
43
44 write(unit=stdout,fmt=*)
45
46 call da_array_print( 1, e, 'Global Eigenvectors' )
47
48 !-------------------------------------------------------------------------
49 ! [2] Test orthogonality of eigenvectors - sum_k (e_m(k) e_n(k)) = delta_mn:
50 !-------------------------------------------------------------------------
51
52 write(unit=stdout,fmt='(A)')' Eigenvector orthogonality check:'
53 write(unit=stdout,fmt='(A)')' Mode Diagonal Maximum off-diagonal'
54
55 do k1 = 1, kz
56 do k2 = 1, kz
57 work(k1,k2) = sum(e(1:kz,k1) * e(1:kz,k2))
58 end do
59
60 array_mask(1:kz) =.true.
61 array_mask(k1) = .false.
62 max_off_diag = maxval(abs(work(k1,:)),mask=array_mask(:))
63 write(unit=stdout,fmt='(I4,4x,1pe12.4,10x,1pe12.4)')k1, work(k1,k1), max_off_diag
64 end do
65 write(unit=stdout,fmt=*)
66
67 !-------------------------------------------------------------------------
68 ! [3] Test eigenvectors completeness - sum_m (e_m(k1) e_m(k2)) = delta_k1k2:
69 !-------------------------------------------------------------------------
70
71 write(unit=stdout,fmt='(A)')' Eigenvector completeness check:'
72 write(unit=stdout,fmt='(A)')' Level Diagonal Maximum off-diagonal'
73
74 do k1 = 1, kz
75 do k2 = 1, kz
76 work(k1,k2) = sum(e(k1,1:kz) * e(k2,1:kz))
77 end do
78
79 array_mask(1:kz) =.true.
80 array_mask(k1) = .false.
81 max_off_diag = maxval(abs(work(k1,:)),mask=array_mask(:))
82 write(unit=stdout,fmt='(I4,4x,1pe12.4,10x,1pe12.4)')k1, work(k1,k1), max_off_diag
83 end do
84 write(unit=stdout,fmt=*)
85
86 !-------------------------------------------------------------------------
87 ! [4] check B correctness: B = E*L*E^T
88 !-------------------------------------------------------------------------
89
90 write(unit=stdout,fmt='(a/a)') &
91 'real and Calculated B (diagonal)', &
92 'lvl real-B Calculated-B'
93
94 do k=1,kz
95 do m=1,kz
96 work(k,m)=l(k)*e(m,k)
97 bc(k,m)=0.0
98 end do
99 end do
100
101 do k1=1,kz
102 do k2=1,kz
103 do m=1,kz
104 bc(k1,k2)=bc(k1,k2)+e(k1,m)*work(m,k2)
105 end do
106 end do
107
108 write(unit=stdout,fmt='(I5,2F20.5)') k1, bx(k1,k1), bc(k1,k1)
109 end do
110
111 do k2=1,kz
112 write(unit=stdout, fmt='(a,i4/a)') &
113 'real and Calculated B (off diagonal):', k2, &
114 'lvl real-B Calculated-B'
115
116 do k1=1,kz
117 write(unit=stdout,fmt='(I5,2F20.5)') k1, bx(k1,k2), bc(k1,k2)
118 end do
119 end do
120
121 if (trace_use) call da_trace_exit("da_eof_decomposition_test")
122
123 end subroutine da_eof_decomposition_test
124
125