FMS 2025.01.02-dev
Flexible Modeling System
Loading...
Searching...
No Matches
fms_diag_reduction_methods.inc
1!***********************************************************************
2!* GNU Lesser General Public License
3!*
4!* This file is part of the GFDL Flexible Modeling System (FMS).
5!*
6!* FMS is free software: you can redistribute it and/or modify it under
7!* the terms of the GNU Lesser General Public License as published by
8!* the Free Software Foundation, either version 3 of the License, or (at
9!* your option) any later version.
10!*
11!* FMS is distributed in the hope that it will be useful, but WITHOUT
12!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14!* for more details.
15!*
16!* You should have received a copy of the GNU Lesser General Public
17!* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
18!***********************************************************************
19
20! for any debug prints
21#ifndef DEBUG_REDUCT
22#define DEBUG_REDUCT .false.
23#endif
24
25!> @brief Do the time_none reduction method (i.e copy the correct portion of the input data)
26subroutine do_time_none_ (data_out, data_in, mask, is_masked, bounds_in, bounds_out, missing_value)
27 real(FMS_TRM_KIND_), intent(inout) :: data_out(:,:,:,:,:) !< output data
28 real(FMS_TRM_KIND_), intent(in) :: data_in(:,:,:,:) !< data to update the buffer with
29 logical, intent(in) :: mask(:,:,:,:) !< mask
30 logical, intent(in) :: is_masked !< .True. if the field is using a mask
31 type(fmsDiagIbounds_type), intent(in) :: bounds_in !< indices indicating the correct portion
32 !! of the input buffer
33 type(fmsDiagIbounds_type), intent(in) :: bounds_out !< indices indicating the correct portion
34 !! of the output buffer
35 real(FMS_TRM_KIND_), intent(in) :: missing_value !< Missing_value for data points that are masked
36
37 integer :: is_in, ie_in, js_in, je_in, ks_in, ke_in !< Starting and ending indices of each dimention for
38 !! the input buffer
39 integer :: is_out, ie_out, js_out, je_out, ks_out, ke_out !< Starting and ending indices of each dimention for
40 !! the output buffer
41
42 is_out = bounds_out%get_imin()
43 ie_out = bounds_out%get_imax()
44 js_out = bounds_out%get_jmin()
45 je_out = bounds_out%get_jmax()
46 ks_out = bounds_out%get_kmin()
47 ke_out = bounds_out%get_kmax()
48
49 is_in = bounds_in%get_imin()
50 ie_in = bounds_in%get_imax()
51 js_in = bounds_in%get_jmin()
52 je_in = bounds_in%get_jmax()
53 ks_in = bounds_in%get_kmin()
54 ke_in = bounds_in%get_kmax()
55
56 if (is_masked) then
57 where (mask(is_in:ie_in, js_in:je_in, ks_in:ke_in, :))
58 data_out(is_out:ie_out, js_out:je_out, ks_out:ke_out, :, 1) = &
59 data_in(is_in:ie_in, js_in:je_in, ks_in:ke_in, :)
60 elsewhere
61 data_out(is_out:ie_out, js_out:je_out, ks_out:ke_out, :, 1) = missing_value
62 end where
63 else
64 data_out(is_out:ie_out, js_out:je_out, ks_out:ke_out, :, 1) = &
65 data_in(is_in:ie_in, js_in:je_in, ks_in:ke_in, :)
66 endif
67
68end subroutine do_time_none_
69
70!> @brief Do the time_min reduction method (i.e maintain the minimum value of the averaging time)
71subroutine do_time_min_ (data_out, data_in, mask, is_masked, bounds_in, bounds_out, missing_value)
72 real(FMS_TRM_KIND_), intent(inout) :: data_out(:,:,:,:,:) !< output data
73 real(FMS_TRM_KIND_), intent(in) :: data_in(:,:,:,:) !< data to update the buffer with
74 logical, intent(in) :: mask(:,:,:,:) !< mask
75 logical, intent(in) :: is_masked !< .True. if the field is using a mask
76 type(fmsDiagIbounds_type), intent(in) :: bounds_in !< indices indicating the correct portion
77 !! of the input buffer
78 type(fmsDiagIbounds_type), intent(in) :: bounds_out !< indices indicating the correct portion
79 !! of the output buffer
80 real(FMS_TRM_KIND_), intent(in) :: missing_value !< Missing_value for data points that are masked
81
82 integer :: is_in, ie_in, js_in, je_in, ks_in, ke_in !< Starting and ending indices of each dimention for
83 !! the input buffer
84 integer :: is_out, ie_out, js_out, je_out, ks_out, ke_out !< Starting and ending indices of each dimention for
85 !! the output buffer
86
87 integer :: i, j, k, l !< For looping
88
89 is_out = bounds_out%get_imin()
90 ie_out = bounds_out%get_imax()
91 js_out = bounds_out%get_jmin()
92 je_out = bounds_out%get_jmax()
93 ks_out = bounds_out%get_kmin()
94 ke_out = bounds_out%get_kmax()
95
96 is_in = bounds_in%get_imin()
97 ie_in = bounds_in%get_imax()
98 js_in = bounds_in%get_jmin()
99 je_in = bounds_in%get_jmax()
100 ks_in = bounds_in%get_kmin()
101 ke_in = bounds_in%get_kmax()
102
103 !> Seperated this loops for performance. If is_masked = .false. (i.e "mask" and "rmask" were never passed in)
104 !! then mask will always be .True. so the if (mask) is redudant.
105 if (is_masked) then
106 do l = 0, size(data_out, 4) - 1
107 do k = 0, ke_out - ks_out
108 do j = 0, je_out - js_out
109 do i = 0, ie_out - is_out
110 if (mask(is_in + i, js_in + j, ks_in + k, l + 1)) then
111 if (data_out(is_out + i, js_out + j, ks_out + k, l + 1, 1) .gt. &
112 data_in(is_in + i, js_in + j, ks_in + k, l + 1) ) then
113 data_out(is_out + i, js_out + j, ks_out + k, l + 1, 1) = &
114 data_in(is_in +i, js_in + j, ks_in + k, l + 1)
115 endif
116 else
117 data_out(is_out + i, js_out + j, ks_out + k, l + 1, 1) = missing_value
118 endif
119 enddo
120 enddo
121 enddo
122 enddo
123 else
124 do l = 0, size(data_out, 4) - 1
125 do k = 0, ke_out - ks_out
126 do j = 0, je_out - js_out
127 do i = 0, ie_out - is_out
128 if (data_out(is_out + i, js_out + j, ks_out + k, l + 1, 1) .gt. &
129 data_in(is_in + i, js_in + j, ks_in + k, l + 1) ) then
130 data_out(is_out + i, js_out + j, ks_out + k, l + 1, 1) = &
131 data_in(is_in +i, js_in + j, ks_in + k, l + 1)
132 endif
133 enddo
134 enddo
135 enddo
136 enddo
137 endif
138
139end subroutine do_time_min_
140
141!> @brief Do the time_max reduction method (i.e maintain the maximum value of the averaging time)
142subroutine do_time_max_ (data_out, data_in, mask, is_masked, bounds_in, bounds_out, missing_value)
143 real(FMS_TRM_KIND_), intent(inout) :: data_out(:,:,:,:,:) !< output data
144 real(FMS_TRM_KIND_), intent(in) :: data_in(:,:,:,:) !< data to update the buffer with
145 logical, intent(in) :: mask(:,:,:,:) !< mask
146 logical, intent(in) :: is_masked !< .True. if the field is using a mask
147 type(fmsDiagIbounds_type), intent(in) :: bounds_in !< indices indicating the correct portion
148 !! of the input buffer
149 type(fmsDiagIbounds_type), intent(in) :: bounds_out !< indices indicating the correct portion
150 !! of the output buffer
151 real(FMS_TRM_KIND_), intent(in) :: missing_value !< Missing_value for data points that are masked
152
153 integer :: is_in, ie_in, js_in, je_in, ks_in, ke_in !< Starting and ending indices of each dimention for
154 !! the input buffer
155 integer :: is_out, ie_out, js_out, je_out, ks_out, ke_out !< Starting and ending indices of each dimention for
156 !! the output buffer
157
158 integer :: i, j, k, l !< For looping
159
160 is_out = bounds_out%get_imin()
161 ie_out = bounds_out%get_imax()
162 js_out = bounds_out%get_jmin()
163 je_out = bounds_out%get_jmax()
164 ks_out = bounds_out%get_kmin()
165 ke_out = bounds_out%get_kmax()
166
167 is_in = bounds_in%get_imin()
168 ie_in = bounds_in%get_imax()
169 js_in = bounds_in%get_jmin()
170 je_in = bounds_in%get_jmax()
171 ks_in = bounds_in%get_kmin()
172 ke_in = bounds_in%get_kmax()
173
174 !> Seperated this loops for performance. If is_masked = .false. (i.e "mask" and "rmask" were never passed in)
175 !! then mask will always be .True. so the if (mask) is redudant.
176 if (is_masked) then
177 do l = 0, size(data_out, 4) - 1
178 do k = 0, ke_out - ks_out
179 do j = 0, je_out - js_out
180 do i = 0, ie_out - is_out
181 if (mask(is_in + i, js_in + j, ks_in + k, l + 1)) then
182 if (data_out(is_out + i, js_out + j, ks_out + k, l + 1, 1) .lt. &
183 data_in(is_in + i, js_in + j, ks_in + k, l + 1) ) then
184 data_out(is_out + i, js_out + j, ks_out + k, l + 1, 1) = &
185 data_in(is_in +i, js_in + j, ks_in + k, l + 1)
186 endif
187 else
188 data_out(is_out + i, js_out + j, ks_out + k, l + 1, 1) = missing_value
189 endif
190 enddo
191 enddo
192 enddo
193 enddo
194 else
195 do l = 0, size(data_out, 4) - 1
196 do k = 0, ke_out - ks_out
197 do j = 0, je_out - js_out
198 do i = 0, ie_out - is_out
199 if (data_out(is_out + i, js_out + j, ks_out + k, l + 1, 1) .lt. &
200 data_in(is_in + i, js_in + j, ks_in + k, l + 1) ) then
201 data_out(is_out + i, js_out + j, ks_out + k, l + 1, 1) = &
202 data_in(is_in +i, js_in + j, ks_in + k, l + 1)
203 endif
204 enddo
205 enddo
206 enddo
207 enddo
208 endif
209end subroutine do_time_max_
210
211!> Update the output buffer for reductions that involve summation (sum, avg, rms, pow).
212!! Elements of the running field output buffer (data_out) are set with the following:
213!!
214!! buffer(l) = buffer(l) + (weight * field(l)) ^ pow
215!!
216!! Where l are the indices passed in through the bounds_in/out
217subroutine do_time_sum_update_(data_out, weight_sum, data_in, mask, is_masked, mask_variant, bounds_in, bounds_out, &
218 missing_value, diurnal_section, weight, pow)
219 real(FMS_TRM_KIND_), intent(inout) :: data_out(:,:,:,:,:) !< output data
220 real(r8_kind), intent(inout) :: weight_sum(:,:,:,:) !< Sum of weights from the output buffer object
221 real(FMS_TRM_KIND_), intent(in) :: data_in(:,:,:,:) !< data to update the buffer with
222 logical, intent(in) :: mask(:,:,:,:) !< mask
223 logical, intent(in) :: is_masked !< .True. if the field is using a mask
224 logical, intent(in) :: mask_variant !< .True. if the mask changes over time
225 type(fmsDiagIbounds_type), intent(in) :: bounds_in !< indices indicating the correct portion
226 !! of the input buffer
227 type(fmsDiagIbounds_type), intent(in) :: bounds_out !< indices indicating the correct portion
228 !! of the output buffer
229 real(FMS_TRM_KIND_), intent(in) :: missing_value !< Missing_value for data points that are masked
230 integer, intent(in) :: diurnal_section !< the diurnal "section" if doing a diurnal reduction
231 !! indicates which index to add data on 5th axis
232 !! if not doing a diurnal reduction, this should always =1
233 real(r8_kind),optional, intent(in) :: weight !< Weight applied to data_in before added to data_out
234 !! used for weighted averages, default 1.0
235 integer ,optional, intent(in) :: pow !< Used for pow(er) reduction,
236 !! calculates field_data^pow before adding to buffer
237
238 real(FMS_TRM_KIND_) :: weight_scale !< local copy of optional weight
239 integer, parameter :: kindl = fms_trm_kind_ !< real kind size as set by macro
240 integer :: diurnal !< diurnal index to indicate which daily section is updated
241 !! will be 1 unless using a diurnal reduction
242
243 if(present(weight)) then
244 weight_scale = real(weight, kind=kindl)
245 else
246 weight_scale = 1.0_kindl
247 endif
248
249 if(diurnal_section .lt. 0) then
250 diurnal = 1
251 else
252 diurnal = diurnal_section
253 endif
254
255 if (is_masked) then
256 if (mask_variant) then
257 ! Mask changes over time so the weight is an array
258 call sum_mask_variant(data_out, data_in, weight_sum, bounds_in, bounds_out, mask, diurnal, weight_scale, pow)
259 else
260 call sum_mask(data_out, data_in, weight_sum, bounds_in, bounds_out, mask, diurnal, &
261 missing_value, weight_scale, pow)
262 endif
263 else
264 call sum_no_mask(data_out, data_in, weight_sum, bounds_in, bounds_out, diurnal, weight_scale, pow)
265 endif
266end subroutine do_time_sum_update_
267
268subroutine sum_mask_(data_out, data_in, weight_sum, bounds_in, bounds_out, mask, diurnal, missing_value, &
269 weight_scale, pow)
270 real(FMS_TRM_KIND_), intent(inout) :: data_out(:,:,:,:,:) !< output data
271 real(FMS_TRM_KIND_), intent(in) :: data_in(:,:,:,:) !< data to update the buffer with
272 real(r8_kind), intent(inout) :: weight_sum(:,:,:,:) !< Sum of weights from the output buffer object
273 type(fmsDiagIbounds_type), intent(in) :: bounds_in !< indices indicating the correct portion
274 !! of the input buffer
275 type(fmsDiagIbounds_type), intent(in) :: bounds_out !< indices indicating the correct portion
276 !! of the output buffer
277 logical, intent(in) :: mask(:,:,:,:) !< mask
278 integer, intent(in) :: diurnal !< diurnal index to indicate which daily section is
279 !! updated will be 1 unless using a diurnal reduction
280 real(FMS_TRM_KIND_), intent(in) :: missing_value !< Missing_value for data points that are masked
281 real(FMS_TRM_KIND_), intent(in) :: weight_scale !< weight scale to use
282 integer ,optional, intent(in) :: pow !< Used for pow(er) reduction,
283 !! calculates field_data^pow before adding to buffer
284
285 integer :: is_in, ie_in, js_in, je_in, ks_in, ke_in !< Starting and ending indices of each dimention for
286 !! the input buffer
287 integer :: is_out, ie_out, js_out, je_out, ks_out, ke_out !< Starting and ending indices of each dimention for
288 !! the output buffer
289 integer :: pow_loc !> local copy of optional pow value (set if using pow reduction)
290 integer :: i, j, k, l !< For looping
291
292 is_out = bounds_out%get_imin()
293 ie_out = bounds_out%get_imax()
294 js_out = bounds_out%get_jmin()
295 je_out = bounds_out%get_jmax()
296 ks_out = bounds_out%get_kmin()
297 ke_out = bounds_out%get_kmax()
298
299 is_in = bounds_in%get_imin()
300 ie_in = bounds_in%get_imax()
301 js_in = bounds_in%get_jmin()
302 je_in = bounds_in%get_jmax()
303 ks_in = bounds_in%get_kmin()
304 ke_in = bounds_in%get_kmax()
305
306 weight_sum = weight_sum + weight_scale
307 if (present(pow)) then
308 do k = 0, ke_out - ks_out
309 do j = 0, je_out - js_out
310 do i = 0, ie_out - is_out
311 where (mask(is_in + i, js_in + j, ks_in + k, :))
312 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) = &
313 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) &
314 + (data_in(is_in +i, js_in + j, ks_in + k, :) * weight_scale) ** pow
315 elsewhere
316 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) = missing_value
317 endwhere
318 enddo
319 enddo
320 enddo
321 else
322 do k = 0, ke_out - ks_out
323 do j = 0, je_out - js_out
324 do i = 0, ie_out - is_out
325 where (mask(is_in + i, js_in + j, ks_in + k, :))
326 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) = &
327 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) &
328 + (data_in(is_in +i, js_in + j, ks_in + k, :) * weight_scale)
329 elsewhere
330 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) = missing_value
331 endwhere
332 enddo
333 enddo
334 enddo
335 endif
336end subroutine sum_mask_
337
338subroutine sum_mask_variant_(data_out, data_in, weight_sum, bounds_in, bounds_out, mask, diurnal, weight_scale, pow)
339 real(FMS_TRM_KIND_), intent(inout) :: data_out(:,:,:,:,:) !< output data
340 real(FMS_TRM_KIND_), intent(in) :: data_in(:,:,:,:) !< data to update the buffer with
341 real(r8_kind), intent(inout) :: weight_sum(:,:,:,:) !< Sum of weights from the output buffer object
342 type(fmsDiagIbounds_type), intent(in) :: bounds_in !< indices indicating the correct portion
343 !! of the input buffer
344 type(fmsDiagIbounds_type), intent(in) :: bounds_out !< indices indicating the correct portion
345 !! of the output buffer
346 logical, intent(in) :: mask(:,:,:,:) !< mask
347 integer, intent(in) :: diurnal !< diurnal index to indicate which daily section is
348 !! updated will be 1 unless using a diurnal reduction
349 real(FMS_TRM_KIND_), intent(in) :: weight_scale !< weight scale to use
350 integer ,optional, intent(in) :: pow !< Used for pow(er) reduction,
351 !! calculates field_data^pow before adding to buffer
352
353 integer :: is_in, ie_in, js_in, je_in, ks_in, ke_in !< Starting and ending indices of each dimention for
354 !! the input buffer
355 integer :: is_out, ie_out, js_out, je_out, ks_out, ke_out !< Starting and ending indices of each dimention for
356 !! the output buffer
357 integer :: pow_loc !> local copy of optional pow value (set if using pow reduction)
358 integer :: i, j, k, l !< For looping
359
360 is_out = bounds_out%get_imin()
361 ie_out = bounds_out%get_imax()
362 js_out = bounds_out%get_jmin()
363 je_out = bounds_out%get_jmax()
364 ks_out = bounds_out%get_kmin()
365 ke_out = bounds_out%get_kmax()
366
367 is_in = bounds_in%get_imin()
368 ie_in = bounds_in%get_imax()
369 js_in = bounds_in%get_jmin()
370 je_in = bounds_in%get_jmax()
371 ks_in = bounds_in%get_kmin()
372 ke_in = bounds_in%get_kmax()
373
374 if (present(pow)) then
375 do k = 0, ke_out - ks_out
376 do j = 0, je_out - js_out
377 do i = 0, ie_out - is_out
378 where (mask(is_in + i, js_in + j, ks_in + k, :))
379 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) = &
380 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) &
381 + (data_in(is_in +i, js_in + j, ks_in + k, :) * weight_scale) ** pow
382
383 !Increase the weight sum for the grid point that was not masked
384 weight_sum(is_out + i, js_out + j, ks_out + k, :) = &
385 weight_sum(is_out + i, js_out + j, ks_out + k, :) + weight_scale
386 endwhere
387 enddo
388 enddo
389 enddo
390 else
391 do k = 0, ke_out - ks_out
392 do j = 0, je_out - js_out
393 do i = 0, ie_out - is_out
394 where (mask(is_in + i, js_in + j, ks_in + k, :))
395 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) = &
396 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) &
397 + (data_in(is_in +i, js_in + j, ks_in + k, :) * weight_scale)
398
399 !Increase the weight sum for the grid point that was not masked
400 weight_sum(is_out + i, js_out + j, ks_out + k, :) = &
401 weight_sum(is_out + i, js_out + j, ks_out + k, :) + weight_scale
402 endwhere
403 enddo
404 enddo
405 enddo
406 endif
407end subroutine sum_mask_variant_
408
409subroutine sum_no_mask_(data_out, data_in, weight_sum, bounds_in, bounds_out, diurnal, weight_scale, pow)
410 real(FMS_TRM_KIND_), intent(inout) :: data_out(:,:,:,:,:) !< output data
411 real(FMS_TRM_KIND_), intent(in) :: data_in(:,:,:,:) !< data to update the buffer with
412 real(r8_kind), intent(inout) :: weight_sum(:,:,:,:) !< Sum of weights from the output buffer object
413 type(fmsDiagIbounds_type), intent(in) :: bounds_in !< indices indicating the correct portion
414 !! of the input buffer
415 type(fmsDiagIbounds_type), intent(in) :: bounds_out !< indices indicating the correct portion
416 !! of the output buffer
417 integer, intent(in) :: diurnal !< diurnal index to indicate which daily section is
418 !! updated will be 1 unless using a diurnal reduction
419 real(FMS_TRM_KIND_), intent(in) :: weight_scale !< weight scale to use
420 integer ,optional, intent(in) :: pow !< Used for pow(er) reduction,
421 !! calculates field_data^pow before adding to buffer
422
423 integer :: is_in, ie_in, js_in, je_in, ks_in, ke_in !< Starting and ending indices of each dimention for
424 !! the input buffer
425 integer :: is_out, ie_out, js_out, je_out, ks_out, ke_out !< Starting and ending indices of each dimention for
426 !! the output buffer
427 integer :: i, j, k, l !< For looping
428
429 is_out = bounds_out%get_imin()
430 ie_out = bounds_out%get_imax()
431 js_out = bounds_out%get_jmin()
432 je_out = bounds_out%get_jmax()
433 ks_out = bounds_out%get_kmin()
434 ke_out = bounds_out%get_kmax()
435
436 is_in = bounds_in%get_imin()
437 ie_in = bounds_in%get_imax()
438 js_in = bounds_in%get_jmin()
439 je_in = bounds_in%get_jmax()
440 ks_in = bounds_in%get_kmin()
441 ke_in = bounds_in%get_kmax()
442
443 weight_sum = weight_sum + weight_scale
444
445 if (present(pow)) then
446 do k = 0, ke_out - ks_out
447 do j = 0, je_out - js_out
448 do i = 0, ie_out - is_out
449 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) = &
450 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) &
451 + (data_in(is_in +i, js_in + j, ks_in + k, :) * weight_scale) ** pow
452 enddo
453 enddo
454 enddo
455 else
456 do k = 0, ke_out - ks_out
457 do j = 0, je_out - js_out
458 do i = 0, ie_out - is_out
459 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) = &
460 data_out(is_out + i, js_out + j, ks_out + k, :, diurnal) &
461 + (data_in(is_in +i, js_in + j, ks_in + k, :) * weight_scale)
462 enddo
463 enddo
464 enddo
465 endif
466end subroutine sum_no_mask_
467
468!> To be called with diag_send_complete, finishes reductions
469!! Just divides the buffer by the counter array(which is just the sum of the weights used in the buffer's reduction)
470!! TODO: change has_mask to an actual logical mask so we don't have to check for missing values
471subroutine sum_update_done_(out_buffer_data, weight_sum, reduction_method, missing_val, has_mask, mask_variant, &
472 n_diurnal_samples)
473 real(FMS_TRM_KIND_), intent(inout) :: out_buffer_data(:,:,:,:,:) !< data buffer previously updated with
474 !! do_time_sum_update
475 real(r8_kind), intent(in) :: weight_sum(:,:,:,:) !< sum of weights for averaging,
476 !! provided via argument to send data
477 integer, intent(in) :: reduction_method !< which reduction method to use
478 !! should always be one of time_avg, time_diurnal, or time_rms
479 real(FMS_TRM_KIND_), intent(in) :: missing_val !< missing value for masked elements
480 logical, intent(in) :: has_mask !< indicates if mask is used so missing values can be skipped
481 logical, intent(in) :: mask_variant !< Indicates if the mask changes over time
482 integer, optional, intent(in) :: n_diurnal_samples !< number of diurnal samples as set in reduction method
483 integer, allocatable :: wsum(:,:,:,:) !< local cp of weight_sum, only changed if using diurnal
484 !! TODO replace conditional in the `where` with passed in and ajusted mask from the original call
485 !logical, optional, intent(in) :: mask(:,:,:,:) !< logical mask from accept data call, if using one.
486 !logical :: has_mask !< whether or not mask is present
487
488 integer :: i, j, k, l !< For do loops
489
490 allocate(wsum(size(weight_sum,1), size(weight_sum,3), size(weight_sum,3), size(weight_sum,4)))
491 ! need to divide weight sum by amount of samples to get the actual
492 ! number of times that the diurnal section was incremented
493 ! legacy diag manager stored these weights explicitly, this doesn't so assumes uniformity in when data is sent
494 if(reduction_method .eq. time_diurnal) then
495 if(.not. present(n_diurnal_samples)) call mpp_error(fatal, &
496 "SUM_UPDATE_DONE_ :: reduction method is diurnal but no sample size was given")
497 wsum = weight_sum / n_diurnal_samples
498 else
499 wsum = weight_sum
500 endif
501
502 if ( has_mask ) then
503 if (.not. mask_variant) then
504 ! The mask does not change over time so wsum is just an integer and it is the same value for all fields
505 where(out_buffer_data(:,:,:,:,:) .ne. missing_val)
506 out_buffer_data(:,:,:,:,:) = out_buffer_data(:,:,:,:,:) &
507 / wsum(1,1,1,1)
508 endwhere
509 else
510 ! The mask changes over time
511 do l = 1, size(out_buffer_data, 4)
512 do k = 1, size(out_buffer_data, 3)
513 do j = 1, size(out_buffer_data, 2)
514 do i = 1, size(out_buffer_data, 1)
515 if (wsum(i, j, k, l) .gt. 0) then
516 out_buffer_data(i,j,k,l,:) = out_buffer_data(i,j,k,l,:)/ wsum(i,j,k,l)
517 else
518 ! Data was never received
519 out_buffer_data(i,j,k,l,:) = missing_val
520 endif
521 enddo
522 enddo
523 enddo
524 enddo
525 endif
526 else
527 ! There is no mask!
528 out_buffer_data(:,:,:,:,:) = out_buffer_data(:,:,:,:,:) &
529 / wsum(1,1,1,1)
530 endif
531
532 if(reduction_method .eq. time_rms .and. has_mask) then
533 where(out_buffer_data(:,:,:,:,1) .ne. missing_val)
534 out_buffer_data(:,:,:,:,1) = sqrt(out_buffer_data(:,:,:,:,1))
535 endwhere
536 else if(reduction_method .eq. time_rms) then
537 out_buffer_data(:,:,:,:,1) = sqrt(out_buffer_data(:,:,:,:,1))
538 endif
539
540end subroutine
541