Line data Source code
1 : !> \file mo_corr.f90
2 : !> \brief \copybrief mo_corr
3 : !> \details \copydetails mo_corr
4 :
5 : !> \brief Provides autocorrelation function calculations.
6 : !> \author Sebastian Mueller
7 : !> \date Dec 2019
8 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
9 : !! FORCES is released under the LGPLv3+ license \license_note
10 : MODULE mo_corr
11 :
12 : USE mo_kind, ONLY : i4, sp, dp
13 :
14 : Implicit NONE
15 :
16 : PUBLIC :: autocorr ! Autocorrelation coefficient at lag k = autocoeffk(k)/autocoeffk(0)
17 :
18 : ! ------------------------------------------------------------------
19 :
20 : !> \brief Autocorrelation function with lag k.
21 :
22 : !> \details Element at lag k of autocorrelation function
23 : !! \f[ ACF(x,k) = \frac{s_{x,k}}{s_{x,0}} \f]
24 : !! where \f$ s_{x,k} \f$ is the autocorrelation coefficient
25 : !! \f[ s_{x,k} = \langle (x_i-\bar{x})(x_{i+k}-\bar{x})\rangle \f]
26 : !!
27 : !! If an optional mask is given, the calculations are only over those locations that correspond
28 : !! to true values in the mask.\n
29 : !! \f$ x \f$ can be single or double precision. The result will have the same numerical precision.
30 : !!
31 : !! \b Example
32 : !!
33 : !! Autocorrelation of 0 time steps
34 : !! \code{.f90}
35 : !! ak = autocorr(x, 0, mask=mask)
36 : !! ---> ak = 1
37 : !! \endcode
38 :
39 : !> \param[in] "real(sp/dp) :: x(:)" Time series.
40 : !> \param[in] "integer(i4) :: k[(:)]" Lag for autocorrelation.
41 : !> \param[in] "optional, logical :: mask(:)" 1D-array of logical values with `size(vec)`.
42 : !! If present, only those locations in vec corresponding
43 : !! to the true values in mask are used.
44 : !> \retval "real(sp/dp) :: ak[(:)]" Coefficient of autocorrelation function at lag k.
45 :
46 : !> \author Matthias Cuntz
47 : !> \date Nov 2011
48 :
49 : !> \author Stephan Thober
50 : !> \date Nov 2012
51 : !! - added 1d version
52 :
53 : !> \author Sebastion Mueller
54 : !> \date Dec 2019
55 : !! - rewritten
56 :
57 : INTERFACE autocorr
58 : MODULE PROCEDURE autocorr_sp, autocorr_dp, autocorr_1d_sp, autocorr_1d_dp
59 : END INTERFACE autocorr
60 :
61 : ! ------------------------------------------------------------------
62 :
63 : CONTAINS
64 :
65 : ! ------------------------------------------------------------------
66 :
67 10 : FUNCTION autocorr_dp(x, k, mask) result(acf)
68 :
69 : USE mo_moment, ONLY: moment
70 :
71 : IMPLICIT NONE
72 :
73 : REAL(dp), DIMENSION(:), INTENT(IN) :: x
74 : INTEGER(i4), INTENT(IN) :: k
75 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
76 : REAL(dp) :: acf
77 :
78 : INTEGER(i4) :: kk ! absolute value of lag k
79 : INTEGER(i4) :: nn ! number of true values in mask
80 5 : REAL(dp) :: n ! real of nn
81 5 : REAL(dp) :: mean, var ! moments
82 7505 : REAL(dp), DIMENSION(size(x)) :: x_shift
83 5 : LOGICAL, DIMENSION(size(x)) :: maske
84 :
85 7505 : maske(:) = .true.
86 5 : if (present(mask)) then
87 3 : if (size(mask) /= size(x)) stop 'Error autocorr_dp: size(mask) /= size(x)'
88 4506 : maske = mask
89 : end if
90 :
91 : ! calculate mean and population variance of x
92 5 : call moment(x, mean=mean, variance=var, mask=maske, sample=.false.)
93 : ! shift x to 0 mean
94 7505 : x_shift = x - mean
95 : ! use absolute value of k
96 5 : kk = abs(k)
97 5 : nn = size(x)
98 7495 : n = real(count(maske(1 : nn - kk).and.maske(1 + kk : nn)), dp)
99 : ! calculate the auto-correlation function
100 7495 : acf = sum(x_shift(1 : nn - kk) * x_shift(1 + kk : nn), mask = (maske(1 : nn - kk) .and. maske(1 + kk : nn))) / n / var
101 :
102 5 : END FUNCTION autocorr_dp
103 :
104 10 : FUNCTION autocorr_sp(x, k, mask) result(acf)
105 :
106 5 : USE mo_moment, ONLY: moment
107 :
108 : IMPLICIT NONE
109 :
110 : REAL(sp), DIMENSION(:), INTENT(IN) :: x
111 : INTEGER(i4), INTENT(IN) :: k
112 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
113 : REAL(sp) :: acf
114 :
115 : INTEGER(i4) :: kk ! absolute value of lag k
116 : INTEGER(i4) :: nn ! number of true values in mask
117 5 : REAL(sp) :: n ! real of nn
118 5 : REAL(sp) :: mean, var ! moments
119 7505 : REAL(sp), DIMENSION(size(x)) :: x_shift
120 5 : LOGICAL, DIMENSION(size(x)) :: maske
121 :
122 7505 : maske(:) = .true.
123 5 : if (present(mask)) then
124 3 : if (size(mask) /= size(x)) stop 'Error autocorr_sp: size(mask) /= size(x)'
125 4506 : maske = mask
126 : end if
127 :
128 : ! calculate mean and population variance of x
129 5 : call moment(x, mean=mean, variance=var, mask=maske, sample=.false.)
130 : ! shift x to 0 mean
131 7505 : x_shift = x - mean
132 : ! use absolute value of k
133 5 : kk = abs(k)
134 5 : nn = size(x)
135 7495 : n = real(count(maske(1 : nn - kk).and.maske(1 + kk : nn)), sp)
136 : ! calculate the auto-correlation function
137 7495 : acf = sum(x_shift(1 : nn - kk) * x_shift(1 + kk : nn), mask = (maske(1 : nn - kk) .and. maske(1 + kk : nn))) / n / var
138 :
139 5 : END FUNCTION autocorr_sp
140 :
141 10 : FUNCTION autocorr_1d_dp(x, k, mask) result(acf)
142 :
143 : IMPLICIT NONE
144 :
145 : REAL(dp), DIMENSION(:), INTENT(IN) :: x
146 : INTEGER(i4), DIMENSION(:), INTENT(IN) :: k
147 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
148 : INTEGER(i4) :: i
149 : REAL(dp), DIMENSION(size(k)) :: acf
150 :
151 2 : if (present(mask)) then
152 3 : do i = 1, size(k)
153 3 : acf(i) = autocorr(x, k(i), mask)
154 : end do
155 : else
156 3 : do i = 1, size(k)
157 3 : acf(i) = autocorr(x, k(i))
158 : end do
159 : end if
160 :
161 :
162 7 : END FUNCTION autocorr_1d_dp
163 :
164 8 : FUNCTION autocorr_1d_sp(x, k, mask) result(acf)
165 :
166 : IMPLICIT NONE
167 :
168 : REAL(sp), DIMENSION(:), INTENT(IN) :: x
169 : INTEGER(i4), DIMENSION(:), INTENT(IN) :: k
170 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
171 : INTEGER(i4) :: i
172 : REAL(sp), DIMENSION(size(k)) :: acf
173 :
174 2 : if (present(mask)) then
175 3 : do i = 1, size(k)
176 3 : acf(i) = autocorr(x, k(i), mask)
177 : end do
178 : else
179 3 : do i = 1, size(k)
180 3 : acf(i) = autocorr(x, k(i))
181 : end do
182 : end if
183 :
184 4 : END FUNCTION autocorr_1d_sp
185 :
186 : END MODULE mo_corr
|
