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 !-------------------------------------------------------------------------
28 ! [1] Print eigenvalues:
29 !-------------------------------------------------------------------------
30
31 tot_variance = sum(l(1:kz))
32 cum_variance = 0.0
33
34 write(unit=stdout,fmt='(A)')' Mode Eigenvalue Cumulative Variance e(k,k)'
35
36 do k = 1, kz
37 cum_variance = cum_variance + l(k)
38 write(unit=stdout,fmt='(I4,4x,e12.4,10x,f8.4,4x,e12.4)') &
39 k, l(k), cum_variance / tot_variance, e(k,k)
40 end do
41
42 write(unit=stdout,fmt=*)
43
44 call da_array_print( 1, e, 'Global Eigenvectors' )
45
46 !-------------------------------------------------------------------------
47 ! [2] Test orthogonality of eigenvectors - sum_k (e_m(k) e_n(k)) = delta_mn:
48 !-------------------------------------------------------------------------
49
50 write(unit=stdout,fmt='(A)')' Eigenvector orthogonality check:'
51 write(unit=stdout,fmt='(A)')' Mode Diagonal Maximum off-diagonal'
52
53 do k1 = 1, kz
54 do k2 = 1, kz
55 work(k1,k2) = sum(e(1:kz,k1) * e(1:kz,k2))
56 end do
57
58 array_mask(1:kz) =.true.
59 array_mask(k1) = .false.
60 max_off_diag = maxval(abs(work(k1,:)),mask=array_mask(:))
61 write(unit=stdout,fmt='(I4,4x,1pe12.4,10x,1pe12.4)')k1, work(k1,k1), max_off_diag
62 end do
63 write(unit=stdout,fmt=*)
64
65 !-------------------------------------------------------------------------
66 ! [3] Test eigenvectors completeness - sum_m (e_m(k1) e_m(k2)) = delta_k1k2:
67 !-------------------------------------------------------------------------
68
69 write(unit=stdout,fmt='(A)')' Eigenvector completeness check:'
70 write(unit=stdout,fmt='(A)')' Level Diagonal Maximum off-diagonal'
71
72 do k1 = 1, kz
73 do k2 = 1, kz
74 work(k1,k2) = sum(e(k1,1:kz) * e(k2,1:kz))
75 end do
76
77 array_mask(1:kz) =.true.
78 array_mask(k1) = .false.
79 max_off_diag = maxval(abs(work(k1,:)),mask=array_mask(:))
80 write(unit=stdout,fmt='(I4,4x,1pe12.4,10x,1pe12.4)')k1, work(k1,k1), max_off_diag
81 end do
82 write(unit=stdout,fmt=*)
83
84 !-------------------------------------------------------------------------
85 ! [4] check B correctness: B = E*L*E^T
86 !-------------------------------------------------------------------------
87
88 write(unit=stdout,fmt='(a/a)') &
89 'real and Calculated B (diagonal)', &
90 'lvl real-B Calculated-B'
91
92 do k=1,kz
93 do m=1,kz
94 work(k,m)=l(k)*e(m,k)
95 bc(k,m)=0.0
96 end do
97 end do
98
99 do k1=1,kz
100 do k2=1,kz
101 do m=1,kz
102 bc(k1,k2)=bc(k1,k2)+e(k1,m)*work(m,k2)
103 end do
104 end do
105
106 write(unit=stdout,fmt='(I5,2F20.5)') k1, bx(k1,k1), bc(k1,k1)
107 end do
108
109 do k2=1,kz
110 write(unit=stdout, fmt='(a,i4/a)') &
111 'real and Calculated B (off diagonal):', k2, &
112 'lvl real-B Calculated-B'
113
114 do k1=1,kz
115 write(unit=stdout,fmt='(I5,2F20.5)') k1, bx(k1,k2), bc(k1,k2)
116 end do
117 end do
118
119 end subroutine da_eof_decomposition_test
120
121