-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy patheof_mod.f90
178 lines (142 loc) · 7.2 KB
/
eof_mod.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
!***********************************************************************
! DESCRIPTION:
! Module applies the EOF approach to analysis meteorological field
! of two dimensions
!
! REVISION HISTORY:
! Prototype 01/2015 by Shen Yu, SCC
! Revised 12/2015 by Li Zhenkun, SCC
!
!***********************************************************************
module eof_mod
implicit none
public :: statistic
public :: eof
contains
!----------------------------------------------------------------------
! Public subroutine to compute mean and standard deviation statistic
! of a one-dimension array
!----------------------------------------------------------------------
subroutine statistic(n, array, mean, sdev)
implicit none
!------------------------------------------------------------------
! dummy arguments
!------------------------------------------------------------------
integer, intent(in) :: n
real(4), intent(in) :: array(n)
real(4), intent(out) :: mean, sdev
!------------------------------------------------------------------
! local variables
!------------------------------------------------------------------
integer :: i
mean = 0.
do i = 1, n
mean = mean + array(i)
end do
mean = mean / real(n)
sdev = 0.
do i = 1, n
sdev = sdev + ( array(i) - mean ) ** 2
end do
sdev = sqrt(sdev / real(n))
end subroutine statistic
!----------------------------------------------------------------------
! Public subroutine to do EOF decomposition
!----------------------------------------------------------------------
subroutine eof(grid_num, samp_num, min_dim, trans_opt, undef, field, eigenvalue, eigenvector, pc, ierr)
use eigen_mod
implicit none
!------------------------------------------------------------------
! dummy arguments
!------------------------------------------------------------------
integer, intent(in) :: grid_num, samp_num, min_dim
integer, intent(in) :: trans_opt !trans_opt = 0 for raw, 1 for anomaly, 2 for normalized anomoly
real(4), intent(in) :: undef
real(4), intent(inout) :: field(grid_num, samp_num)
real(4), intent(out) :: eigenvalue(min_dim)
real(4), intent(out) :: eigenvector(grid_num, min_dim)
real(4), intent(out) :: pc(min_dim, samp_num)
integer, intent(out) :: ierr
!------------------------------------------------------------------
! local variables
!------------------------------------------------------------------
character(len=*), parameter :: subname = '(eof)'
integer, parameter :: opt = 3 !opt = 1, 2 find eigenvalues only, opt = 3, 4 find both eigenvalues and its eigenvectors.
real(4) :: symm_matrix(min_dim, min_dim), symm_eigenvector(min_dim, min_dim)
real(4) :: mean(grid_num), sdev(grid_num)
real(4) :: lamda(min_dim)
integer :: abnorm_num
real(4) :: minima = 1.e-25
integer :: i, j
!----------------------------------------------------------------------
! Check whether there are undef values in matrix field
!----------------------------------------------------------------------
ierr = 0
abnorm_num = count(field == undef)
if ( abnorm_num /= 0 ) then
write( 6, '(2a)' ) subname, 'ERROR: Undefind value found in matrix'
ierr = 1
return
end if
!----------------------------------------------------------------------
! Compute mean and standard deviation of matrix field for second
! dimension samp_num
!----------------------------------------------------------------------
do i = 1, grid_num
call statistic(samp_num, field(i, :), mean(i), sdev(i))
end do
!----------------------------------------------------------------------
! Matrix transform according to the given parameter trans_opt
!----------------------------------------------------------------------
if ( trans_opt == 1 ) then ! anomaly transform
forall ( i = 1:samp_num ) field(1:grid_num, i) = field(1:grid_num, i) - mean(1:grid_num)
else if ( trans_opt == 2 ) then ! standardized transform
abnorm_num = count(sdev <= minima) ! checking whether there are zeros in sdev
if ( abnorm_num /= 0 ) then
write( 6, '(2a)' ) subname, 'ERROR: Values very close to Zero found in matrix standrad deviation'
ierr = 1
return
end if
forall ( i = 1:samp_num ) field(1:grid_num, i) = ( field(1:grid_num, i) - mean(1:grid_num) ) / sdev(1:grid_num)
end if
!----------------------------------------------------------------------
! Transpose matrix if grid_num > samp_num (i.e. space dimension >
! time dimension
!----------------------------------------------------------------------
symm_matrix = 0.0
if ( grid_num > samp_num ) then
symm_matrix = matmul(transpose(field), field)
symm_matrix = symm_matrix / grid_num
else
symm_matrix = matmul(field, transpose(field))
symm_matrix = symm_matrix / samp_num
end if
!----------------------------------------------------------------------
! Call subroutine rs() to obtain all eigenvalues & corresponding
! eigenvectors of symmetric matrix
!----------------------------------------------------------------------
call rs(min_dim, min_dim, symm_matrix, eigenvalue, opt, symm_eigenvector, ierr)
if ( ierr /= 0 ) then
write( 6, '(2a, i, a)' ) subname, 'ERROR: can NOT evaluate the ', ierr, ' eigenvalue'
return
end if
!----------------------------------------------------------------------
! Reorder the eigenvalues in a descending order as well as corres-
! ponding eigenvectors
!----------------------------------------------------------------------
eigenvalue(1:min_dim) = eigenvalue(min_dim:1:-1)
symm_eigenvector(:, 1:min_dim) = symm_eigenvector(:, min_dim:1:-1)
if ( grid_num > samp_num ) then ! if grid_num > samp_num find eigenvalue and its eigenvector
eigenvalue = (grid_num / samp_num) * eigenvalue ! eigenvalue
eigenvector = matmul(field, symm_eigenvector) ! eigenvector
forall( i = 1:min_dim ) lamda(i) = sqrt(dot_product(eigenvector(:, i), eigenvector(:, i)))
forall(i = 1:min_dim) eigenvector(:, i) = eigenvector(:, i) / lamda(i) ! normalized treatment.
else
eigenvector = symm_eigenvector
end if
!----------------------------------------------------------------------
! Obtain PCs of matrix field
!----------------------------------------------------------------------
pc = matmul(transpose(eigenvector), field)
end subroutine eof
end module eof_mod