feat: modern diag add time_sum reduction (NOAA-GFDL#1375)

diag_manager/diag_data.F90

@@ -108,7 +108,7 @@ MODULE diag_data_mod
                                            !! to indicate to use the full axis instead of a sub-axis
   INTEGER, PARAMETER :: GLO_REG_VAL_ALT = -1 !< Alternate value used in the region specification of the
                                              !! diag_table to indicate to use the full axis instead of a sub-axis
-  REAL, PARAMETER :: CMOR_MISSING_VALUE = 1.0e20 !< CMOR standard missing value
+  REAL(r8_kind), PARAMETER :: CMOR_MISSING_VALUE = 1.0e20 !< CMOR standard missing value
   INTEGER, PARAMETER :: DIAG_FIELD_NOT_FOUND = -1 !< Return value for a diag_field that isn't found in the diag_table
   INTEGER, PARAMETER :: latlon_gridtype = 1
   INTEGER, PARAMETER :: index_gridtype = 2
@@ -391,9 +391,13 @@ MODULE diag_data_mod
                                          !! the default behavior will do the average between day1 hour3 to day2 hour3
   ! <!-- netCDF variable -->
 
-  REAL :: FILL_VALUE = NF_FILL_REAL !< Fill value used.  Value will be <TT>NF90_FILL_REAL</TT> if using the
+#ifdef use_netCDF
+  REAL(r8_kind) :: FILL_VALUE = NF_FILL_REAL !< Fill value used.  Value will be <TT>NF90_FILL_REAL</TT> if using the
                                     !! netCDF module, otherwise will be 9.9692099683868690e+36.
                                     ! from file /usr/local/include/netcdf.inc
+#else
+  REAL(r8_kind) :: FILL_VALUE = 9.9692099683868690e+36
+#endif
 
   !! @note `pack_size` and `pack_size_str` are set in diag_manager_init depending on how FMS was compiled
   !! if FMS was compiled with default reals as 64bit, it will be set to 1 and "double",
@@ -406,8 +410,8 @@ MODULE diag_data_mod
                                              !! set to "double" or "float"
 
   ! <!-- REAL public variables -->
-  REAL :: EMPTY = 0.0
-  REAL :: MAX_VALUE, MIN_VALUE
+  REAL(r8_kind) :: EMPTY = 0.0
+  REAL(r8_kind) :: MAX_VALUE, MIN_VALUE
 
   ! <!-- Global data for all files -->
   TYPE(time_type) :: diag_init_time !< Time diag_manager_init called.  If init_time not included in

diag_manager/diag_manager.F90

@@ -1443,7 +1443,7 @@ LOGICAL FUNCTION send_data_0d(diag_field_id, field, time, err_msg)
     TYPE(time_type), INTENT(in), OPTIONAL :: time
     CHARACTER(len=*), INTENT(out), OPTIONAL :: err_msg
 
-    REAL :: field_out(1, 1, 1) !< Local copy of field
+    CLASS(*), allocatable :: field_out(:, :, :) !< Local copy of field
 
     ! If diag_field_id is < 0 it means that this field is not registered, simply return
     IF ( diag_field_id <= 0 ) THEN
@@ -1454,9 +1454,23 @@ LOGICAL FUNCTION send_data_0d(diag_field_id, field, time, err_msg)
     ! First copy the data to a three d array with last element 1
     SELECT TYPE (field)
     TYPE IS (real(kind=r4_kind))
-       field_out(1, 1, 1) = field
+       allocate(real(r4_kind) :: field_out(1,1,1))
+       select type(field_out)
+         type is (real(r4_kind))
+           field_out(1, 1, 1) = field
+         class default
+           call error_mesg('diag_manager_mod::send_data_0d', &
+               & 'Error allocating field out as real(r4_kind)', FATAL)
+       end select
     TYPE IS (real(kind=r8_kind))
-       field_out(1, 1, 1) = real(field)
+       allocate(real(r8_kind) :: field_out(1,1,1))
+       select type(field_out)
+         type is (real(r8_kind))
+           field_out(1, 1, 1) = field
+         class default
+           call error_mesg('diag_manager_mod::send_data_0d', &
+               & 'Error allocating field out as real(r8_kind)', FATAL)
+       end select
     CLASS DEFAULT
        CALL error_mesg ('diag_manager_mod::send_data_0d',&
             & 'The field is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
@@ -1476,7 +1490,7 @@ LOGICAL FUNCTION send_data_1d(diag_field_id, field, time, is_in, mask, rmask, ie
     LOGICAL, INTENT(in), DIMENSION(:), OPTIONAL :: mask
     CHARACTER(len=*), INTENT(out), OPTIONAL :: err_msg
 
-    REAL, DIMENSION(SIZE(field(:)), 1, 1) :: field_out !< Local copy of field
+    CLASS(*), ALLOCATABLE :: field_out(:,:,:) !< Local copy of field
     LOGICAL, DIMENSION(SIZE(field(:)), 1, 1) ::  mask_out !< Local copy of mask
 
     ! If diag_field_id is < 0 it means that this field is not registered, simply return
@@ -1486,11 +1500,26 @@ LOGICAL FUNCTION send_data_1d(diag_field_id, field, time, is_in, mask, rmask, ie
     END IF
 
     ! First copy the data to a three d array with last element 1
+    ! type checking done in diag_send_data
     SELECT TYPE (field)
     TYPE IS (real(kind=r4_kind))
-       field_out(:, 1, 1) = field
+       allocate(real(r4_kind) :: field_out(SIZE(field),1,1))
+       select type(field_out)
+         type is (real(r4_kind))
+           field_out(:, 1, 1) = field
+         class default
+           call error_mesg('diag_manager_mod::send_data_1d', &
+               & 'Error allocating field out as real(r4_kind)', FATAL)
+       end select
     TYPE IS (real(kind=r8_kind))
-       field_out(:, 1, 1) = real(field)
+       allocate(real(r8_kind) :: field_out(SIZE(field),1,1))
+       select type(field_out)
+         type is (real(r8_kind))
+           field_out(:, 1, 1) = field
+         class default
+           call error_mesg('diag_manager_mod::send_data_1d', &
+               & 'Error allocating field out as real(r8_kind)', FATAL)
+       end select
     CLASS DEFAULT
        CALL error_mesg ('diag_manager_mod::send_data_1d',&
             & 'The field is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
@@ -1545,7 +1574,7 @@ LOGICAL FUNCTION send_data_2d(diag_field_id, field, time, is_in, js_in, &
     CLASS(*), INTENT(in), DIMENSION(:,:),OPTIONAL :: rmask
     CHARACTER(len=*), INTENT(out), OPTIONAL :: err_msg
 
-    REAL, DIMENSION(SIZE(field,1),SIZE(field,2),1) :: field_out !< Local copy of field
+    CLASS(*), ALLOCATABLE :: field_out(:,:,:) !< Local copy of field
     LOGICAL, DIMENSION(SIZE(field,1),SIZE(field,2),1) ::  mask_out !< Local copy of mask
 
     ! If diag_field_id is < 0 it means that this field is not registered, simply return
@@ -1557,9 +1586,23 @@ LOGICAL FUNCTION send_data_2d(diag_field_id, field, time, is_in, js_in, &
     ! First copy the data to a three d array with last element 1
     SELECT TYPE (field)
     TYPE IS (real(kind=r4_kind))
-       field_out(:, :, 1) = field
+       allocate(real(r4_kind) :: field_out(SIZE(field,1),SIZE(field,2),1))
+       select type(field_out)
+         type is (real(r4_kind))
+           field_out(:, :, 1) = field
+         class default
+           call error_mesg('diag_manager_mod::send_data_2d', &
+               & 'Error allocating field out as real(r4_kind)', FATAL)
+       end select
     TYPE IS (real(kind=r8_kind))
-       field_out(:, :, 1) = real(field)
+       allocate(real(r8_kind) :: field_out(SIZE(field,1),SIZE(field,2),1))
+       select type(field_out)
+         type is (real(r8_kind))
+           field_out(:, :, 1) = field
+         class default
+           call error_mesg('diag_manager_mod::send_data_2d', &
+               & 'Error allocating field out as real(r8_kind)', FATAL)
+       end select
     CLASS DEFAULT
        CALL error_mesg ('diag_manager_mod::send_data_2d',&
             & 'The field is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)

diag_manager/fms_diag_object.F90

@@ -517,7 +517,7 @@ logical function fms_diag_accept_data (this, diag_field_id, field_data, mask, rm
                                                               !! the calculationslater.  \note This is experimental
   character(len=128)                       :: error_string    !< Store error text
   logical                                  :: data_buffer_is_allocated !< .true. if the data buffer is allocated
-  character(len=128)                       :: field_info      !< String holding info about the field to append to the
+  character(len=256)                       :: field_info      !< String holding info about the field to append to the
                                                               !! error message
   logical, allocatable, dimension(:,:,:,:) :: oor_mask        !< Out of range mask
   real(kind=r8_kind)                       :: field_weight    !< Weight to use when averaging (it will be converted
@@ -899,6 +899,11 @@ function fms_diag_do_reduction(this, field_data, diag_field_id, oor_mask, weight
         return
       endif
     case (time_sum)
+      error_msg = buffer_ptr%do_time_sum_wrapper(field_data, oor_mask, field_ptr%get_mask_variant(), &
+        bounds_in, bounds_out, missing_value)
+      if (trim(error_msg) .ne. "") then
+        return
+      endif
     case (time_average)
     case (time_power)
     case (time_rms)

diag_manager/fms_diag_output_buffer.F90

@@ -34,7 +34,7 @@ module fms_diag_output_buffer_mod
 use fms2_io_mod, only: FmsNetcdfFile_t, write_data, FmsNetcdfDomainFile_t, FmsNetcdfUnstructuredDomainFile_t
 use fms_diag_yaml_mod, only: diag_yaml
 use fms_diag_bbox_mod, only: fmsDiagIbounds_type
-use fms_diag_reduction_methods_mod, only: do_time_none, do_time_min, do_time_max
+use fms_diag_reduction_methods_mod, only: do_time_none, do_time_min, do_time_max, do_time_sum_update
 use fms_diag_time_utils_mod, only: diag_time_inc
 
 implicit none
@@ -48,10 +48,8 @@ module fms_diag_output_buffer_mod
   class(*), allocatable :: buffer(:,:,:,:,:)  !< 5D numeric data array
   integer               :: ndim               !< Number of dimensions for each variable
   integer,  allocatable :: buffer_dims(:)     !< holds the size of each dimension in the buffer
-  real(r8_kind), allocatable :: counter(:,:,:,:,:) !< (x,y,z, time-of-day) used in the time averaging functions
+  real(r8_kind)         :: weight_sum !< (x,y,z, time-of-day) used in the time averaging functions
   integer,  allocatable :: num_elements(:)    !< used in time-averaging
-  real(r8_kind), allocatable :: count_0d(:)        !< used in time-averaging along with
-                                              !! counter which is stored in the child types (bufferNd)
   integer,  allocatable :: axis_ids(:)        !< Axis ids for the buffer
   integer               :: field_id           !< The id of the field the buffer belongs to
   integer               :: yaml_id            !< The id of the yaml id the buffer belongs to
@@ -78,6 +76,7 @@ module fms_diag_output_buffer_mod
   procedure :: do_time_none_wrapper
   procedure :: do_time_min_wrapper
   procedure :: do_time_max_wrapper
+  procedure :: do_time_sum_wrapper
 
 end type fmsDiagOutputBuffer_type
 
@@ -124,9 +123,7 @@ subroutine flush_buffer(this)
   this%yaml_id     = diag_null
   if (allocated(this%buffer))       deallocate(this%buffer)
   if (allocated(this%buffer_dims))  deallocate(this%buffer_dims)
-  if (allocated(this%counter))      deallocate(this%counter)
   if (allocated(this%num_elements)) deallocate(this%num_elements)
-  if (allocated(this%count_0d))     deallocate(this%count_0d)
   if (allocated(this%axis_ids))     deallocate(this%axis_ids)
 end subroutine flush_buffer
 
@@ -154,38 +151,22 @@ subroutine allocate_buffer(this, buff_type, ndim, buff_sizes, field_name, diurna
     type is (integer(kind=i4_kind))
       allocate(integer(kind=i4_kind) :: this%buffer(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4),  &
                                                   & buff_sizes(5)))
-      allocate(this%counter(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
-                                                   & buff_sizes(5)))
-      allocate(this%count_0d(n_samples))
-      this%counter = 0.0_r4_kind
-      this%count_0d = 0.0_r4_kind
+      this%weight_sum = 0.0_r4_kind
       this%buffer_type = i4
     type is (integer(kind=i8_kind))
       allocate(integer(kind=i8_kind) :: this%buffer(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
                                                   & buff_sizes(5)))
-      allocate(this%counter(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
-                                                  & buff_sizes(5)))
-      allocate(this%count_0d(n_samples))
-      this%counter = 0.0_r8_kind
-      this%count_0d = 0.0_r8_kind
+      this%weight_sum = 0.0_r8_kind
       this%buffer_type = i8
     type is (real(kind=r4_kind))
       allocate(real(kind=r4_kind) :: this%buffer(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4),  &
                                                & buff_sizes(5)))
-      allocate(this%counter(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
-                                                & buff_sizes(5)))
-      allocate(this%count_0d(n_samples))
-      this%counter = 0.0_r4_kind
-      this%count_0d = 0.0_r4_kind
+      this%weight_sum = 0.0_r4_kind
       this%buffer_type = r4
     type is (real(kind=r8_kind))
       allocate(real(kind=r8_kind) :: this%buffer(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
                                                & buff_sizes(5)))
-      allocate(this%counter(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
-                                                & buff_sizes(5)))
-      allocate(this%count_0d(n_samples))
-      this%counter = 0.0_r8_kind
-      this%count_0d = 0.0_r8_kind
+      this%weight_sum = 0.0_r8_kind
       this%buffer_type = r8
     class default
        call mpp_error("allocate_buffer", &
@@ -194,7 +175,6 @@ subroutine allocate_buffer(this, buff_type, ndim, buff_sizes, field_name, diurna
   end select
   allocate(this%num_elements(n_samples))
   this%num_elements = 0
-  this%count_0d   = 0
   this%done_with_math = .false.
   allocate(this%buffer_dims(5))
   this%buffer_dims(1) = buff_sizes(1)
@@ -571,5 +551,42 @@ function do_time_max_wrapper(this, field_data, mask, is_masked, bounds_in, bound
       end select
   end select
 end function do_time_max_wrapper
+
+!> @brief Does the time_sum reduction method on the buffer object
+!! @return Error message if the math was not successful
+function do_time_sum_wrapper(this, field_data, mask, is_masked, bounds_in, bounds_out, missing_value) &
+  result(err_msg)
+  class(fmsDiagOutputBuffer_type), intent(inout) :: this                !< buffer object to write
+  class(*),                        intent(in)    :: field_data(:,:,:,:) !< Buffer data for current time
+  type(fmsDiagIbounds_type),       intent(in)    :: bounds_in           !< Indicies for the buffer passed in
+  type(fmsDiagIbounds_type),       intent(in)    :: bounds_out          !< Indicies for the output buffer
+  logical,                         intent(in)    :: mask(:,:,:,:)       !< Mask for the field
+  logical,                         intent(in)    :: is_masked           !< .True. if the field has a mask
+  real(kind=r8_kind),              intent(in)    :: missing_value       !< Missing_value for data points that are masked
+  character(len=50) :: err_msg
+
+  !TODO This will be expanded for integers
+  err_msg = ""
+  select type (output_buffer => this%buffer)
+    type is (real(kind=r8_kind))
+      select type (field_data)
+      type is (real(kind=r8_kind))
+        call do_time_sum_update(output_buffer, this%weight_sum, field_data, mask, is_masked, &
+                                bounds_in, bounds_out, missing_value)
+      class default
+        err_msg="do_time_sum_wrapper::the output buffer and the buffer send in are not of the same type (r8_kind)"
+      end select
+    type is (real(kind=r4_kind))
+      select type (field_data)
+      type is (real(kind=r4_kind))
+        call do_time_sum_update(output_buffer, this%weight_sum, field_data, mask, is_masked, bounds_in, bounds_out, &
+          real(missing_value, kind=r4_kind))
+      class default
+        err_msg="do_time_sum_wrapper::the output buffer and the buffer send in are not of the same type (r4_kind)"
+      end select
+    class default
+      err_msg="do_time_sum_wrapper::the output buffer is not a valid type, must be real(r8_kind) or real(r4_kind)"
+  end select
+end function do_time_sum_wrapper
 #endif
 end module fms_diag_output_buffer_mod

diag_manager/fms_diag_reduction_methods.F90

@@ -35,7 +35,7 @@ module fms_diag_reduction_methods_mod
   private
 
   public :: check_indices_order, init_mask, set_weight
-  public :: do_time_none, do_time_min, do_time_max
+  public :: do_time_none, do_time_min, do_time_max, do_time_sum_update
 
   !> @brief Does the time_none reduction method. See include/fms_diag_reduction_methods.inc
   !TODO This needs to be extended to integers
@@ -55,6 +55,13 @@ module fms_diag_reduction_methods_mod
     module procedure do_time_max_r4, do_time_max_r8
   end interface do_time_max
 
+  !> @brief Sum update updates the buffer for any reductions that involve summation
+  !! (ie. time_sum, avg, rms, pow)
+  !!TODO This needs to be extended to integers
+  interface do_time_sum_update
+    module procedure do_time_sum_update_r4, do_time_sum_update_r8
+  end interface
+
   contains
 
   !> @brief Checks improper combinations of is, ie, js, and je.