FMS  2025.02
Flexible Modeling System
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)
26 subroutine 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 
68 end subroutine do_time_none_
69 
70 !> @brief Do the time_min reduction method (i.e maintain the minimum value of the averaging time)
71 subroutine 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 
139 end subroutine do_time_min_
140 
141 !> @brief Do the time_max reduction method (i.e maintain the maximum value of the averaging time)
142 subroutine 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
209 end 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
217 subroutine 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
266 end subroutine do_time_sum_update_
267 
268 subroutine 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
336 end subroutine sum_mask_
337 
338 subroutine 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
407 end subroutine sum_mask_variant_
408 
409 subroutine 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
466 end 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
471 subroutine 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 
540 end subroutine
541