Soil Water Balance (SWB2)
Loading...
Searching...
No Matches
running_grid_stats.F90
Go to the documentation of this file.
1module running_grid_stats
2
3 use constants_and_conversions, only : true, false
4 use exceptions, only : assert
5 use grid
6 use iso_c_binding, only : c_long_long, c_double, c_int, c_bool
7 implicit none
8
9 private
10
11 public :: running_stats_t
12
13 type running_stats_t
14 type (general_grid_t), pointer :: grd_delta => null()
15 type (general_grid_t), pointer :: grd_delta_n => null()
16 type (general_grid_t), pointer :: grd_delta_n2 => null()
17 type (general_grid_t), pointer :: grd_term1 => null()
18 type (general_grid_t), pointer :: grd_m1 => null()
19 type (general_grid_t), pointer :: grd_m2 => null()
20 type (general_grid_t), pointer :: grd_sum => null()
21 integer (c_long_long) :: n = 0
22 real (c_double) :: nodata_value = -9999.
23 ! type (GENERAL_GRID_T), pointer :: grd_M3 => null()
24 ! type (GENERAL_GRID_T), pointer :: grd_M4 => null()
25
26 contains
27
28 procedure :: clear => clear_sub
29 procedure :: initialize => initialize_sub
30 procedure :: push => push_sub
31 procedure :: mean => calc_mean_sub
32 procedure :: sum => calc_sum_sub
33 procedure :: variance => calc_variance_sub
34 procedure :: std_deviation => calc_std_deviation_sub
35
36 end type running_stats_t
37
38 contains
39
40 subroutine initialize_sub(this, NX, NY, X0, Y0, X1, Y1, nodata_value)
41
42 class(running_stats_t), intent(inout) :: this
43 integer (c_int), intent(in) :: NX
44 integer (c_int), intent(in) :: NY
45 real (c_double), intent(in) :: X0
46 real (c_double), intent(in) :: Y0
47 real (c_double), intent(in) :: X1
48 real (c_double), intent(in) :: Y1
49 real (c_double), intent(in) :: nodata_value
50
51 this%nodata_value = nodata_value
52
53 this%grd_delta => grid_create( inx=nx, &
54 iny=ny, &
55 rx0=x0, &
56 ry0=y0, &
57 rx1=x1, &
58 ry1=y1, &
59 idatatype=grid_datatype_double )
60 this%grd_delta_n => grid_create( inx=nx, &
61 iny=ny, &
62 rx0=x0, &
63 ry0=y0, &
64 rx1=x1, &
65 ry1=y1, &
66 idatatype=grid_datatype_double )
67
68 this%grd_delta_n2 => grid_create( inx=nx, &
69 iny=ny, &
70 rx0=x0, &
71 ry0=y0, &
72 rx1=x1, &
73 ry1=y1, &
74 idatatype=grid_datatype_double )
75
76 this%grd_term1 => grid_create( inx=nx, &
77 iny=ny, &
78 rx0=x0, &
79 ry0=y0, &
80 rx1=x1, &
81 ry1=y1, &
82 idatatype=grid_datatype_double )
83
84 this%grd_M1 => grid_create( inx=nx, &
85 iny=ny, &
86 rx0=x0, &
87 ry0=y0, &
88 rx1=x1, &
89 ry1=y1, &
90 idatatype=grid_datatype_double )
91
92 this%grd_M2 => grid_create( inx=nx, &
93 iny=ny, &
94 rx0=x0, &
95 ry0=y0, &
96 rx1=x1, &
97 ry1=y1, &
98 idatatype=grid_datatype_double )
99
100 this%grd_sum => grid_create( inx=nx, &
101 iny=ny, &
102 rx0=x0, &
103 ry0=y0, &
104 rx1=x1, &
105 ry1=y1, &
106 idatatype=grid_datatype_double )
107
108 end subroutine initialize_sub
109
110!------------------------------------------------------------------------------
111
112 subroutine clear_sub(this)
113
114 class(running_stats_t), intent(inout) :: this
115
116 if ( allocated(this%grd_delta%dpData) ) this%grd_delta%dpData = 0.0_c_double
117 if ( allocated(this%grd_delta_n%dpData) ) this%grd_delta_n%dpData = 0.0_c_double
118 if ( allocated(this%grd_delta_n2%dpData) ) this%grd_delta_n2%dpData = 0.0_c_double
119 if ( allocated(this%grd_term1%dpData) ) this%grd_term1%dpData = 0.0_c_double
120 if ( allocated(this%grd_M1%dpData) ) this%grd_M1%dpData = 0.0_c_double
121 if ( allocated(this%grd_M2%dpData) ) this%grd_M2%dpData = 0.0_c_double
122 if ( allocated(this%grd_sum%dpData) ) this%grd_sum%dpData = 0.0_c_double
123 this%n = 0_c_long_long
124
125 end subroutine clear_sub
126
127!------------------------------------------------------------------------------
128
129 subroutine push_sub( this, grd_data, mask )
130
131 class(running_stats_t), intent(inout) :: this
132 real (c_double), intent(in) :: grd_data(:,:)
133 logical (c_bool), intent(in) :: mask(:,:)
134
135 integer (c_long_long) :: n1
136
137 n1 = this%n
138 this%n = this%n + 1
139
140 where ( mask )
141
142 this%grd_sum%dpData = this%grd_sum%dpData + grd_data
143 this%grd_delta%dpData = grd_data - this%grd_M1%dpData
144 this%grd_delta_n%dpData = this%grd_delta%dpData / real(this%n, c_double)
145 this%grd_delta_n2%dpData = this%grd_delta_n%dpData * this%grd_delta_n%dpData
146 this%grd_term1%dpData = this%grd_delta%dpData * this%grd_delta_n%dpData &
147 * real(n1, c_double)
148 this%grd_M1%dpData = this%grd_M1%dpData + this%grd_delta_n%dpData
149 this%grd_M2%dpData = this%grd_M2%dpData + this%grd_term1%dpData
150
151 else where
152
153 this%grd_sum%dpData = -9999._c_double
154 this%grd_M1%dpData = -9999._c_double
155 this%grd_M2%dpData = -9999._c_double
156
157 end where
158
159 end subroutine push_sub
160
161!------------------------------------------------------------------------------
162
163 subroutine calc_mean_sub(this, dpData)
164
165 class(running_stats_t), intent(inout) :: this
166 real (c_double), dimension(:,:), intent(out) :: dpData
167
168 dpdata = this%grd_M1%dpData
169
170 end subroutine calc_mean_sub
171
172!------------------------------------------------------------------------------
173
174 subroutine calc_sum_sub(this, dpData)
175
176 class(running_stats_t), intent(inout) :: this
177 real (c_double), dimension(:,:), intent(out) :: dpData
178
179 dpdata = this%grd_sum%dpData
180
181 end subroutine calc_sum_sub
182
183!------------------------------------------------------------------------------
184
185 subroutine calc_variance_sub(this, dpData)
186
187 class(running_stats_t), intent(inout) :: this
188 real (c_double), dimension(:,:), intent(out) :: dpData
189
190! call assert(this%n > 1, "Not enough data points processed to calculate variance")
191
192 where ((this%grd_M1%dpData > -9999._c_double) .and. this%n > 1)
193 dpdata = this%grd_M1%dpData / (real(this%n, c_double) - 1.0_c_double)
194 else where
195 dpdata = -9999._c_double
196 end where
197
198 end subroutine calc_variance_sub
199
200!------------------------------------------------------------------------------
201
202 subroutine calc_std_deviation_sub(this, dpData)
203
204 class(running_stats_t), intent(inout) :: this
205 real (c_double), dimension(:,:), allocatable, intent(out) :: dpData
206
207! call assert(this%n > 1, "Not enough data points processed to calculate variance")
208
209 if(this%n > 1) then
210 dpdata = sqrt(this%grd_M1%dpData / (real(this%n, c_double) - 1.0_c_double))
211 else
212 dpdata = 0.0_c_double
213 endif
214
215 end subroutine calc_std_deviation_sub
216
217!------------------------------------------------------------------------------
218
219end module running_grid_stats
constants_and_conversions
This module contains physical constants and convenience functions aimed at performing unit conversion...
Definition constants_and_conversions.F90:8
constants_and_conversions::true
logical(c_bool), parameter, public true
Definition constants_and_conversions.F90:38
constants_and_conversions::false
logical(c_bool), parameter, public false
Definition constants_and_conversions.F90:39
grid
Provides support for input and output of gridded ASCII data, as well as for creation and destruction ...
Definition grid.F90:8
grid::grid_datatype_double
integer(c_int), parameter, public grid_datatype_double
Definition grid.F90:27
running_grid_stats::calc_sum_sub
subroutine calc_sum_sub(this, dpdata)
Definition running_grid_stats.F90:175
running_grid_stats::calc_mean_sub
subroutine calc_mean_sub(this, dpdata)
Definition running_grid_stats.F90:164
running_grid_stats::push_sub
subroutine push_sub(this, grd_data, mask)
Definition running_grid_stats.F90:130
running_grid_stats::calc_std_deviation_sub
subroutine calc_std_deviation_sub(this, dpdata)
Definition running_grid_stats.F90:203
running_grid_stats::calc_variance_sub
subroutine calc_variance_sub(this, dpdata)
Definition running_grid_stats.F90:186
running_grid_stats::clear_sub
subroutine clear_sub(this)
Definition running_grid_stats.F90:113
running_grid_stats::initialize_sub
subroutine initialize_sub(this, nx, ny, x0, y0, x1, y1, nodata_value)
Definition running_grid_stats.F90:41