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