feat: modern diag add time_average reduction method and test (NOAA-GF…

…DL#1421)

diag_manager/diag_data.F90

@@ -114,13 +114,13 @@ MODULE diag_data_mod
   INTEGER, PARAMETER :: index_gridtype = 2
   INTEGER, PARAMETER :: null_gridtype = DIAG_NULL
   INTEGER, PARAMETER :: time_none    = 0 !< There is no reduction method
-  INTEGER, PARAMETER :: time_average = 1 !< The reduction method is avera
-  INTEGER, PARAMETER :: time_rms     = 2 !< The reduction method is rms
-  INTEGER, PARAMETER :: time_max     = 3 !< The reduction method is max
-  INTEGER, PARAMETER :: time_min     = 4 !< The reduction method is min
-  INTEGER, PARAMETER :: time_sum     = 5 !< The reudction method is sum
+  INTEGER, PARAMETER :: time_min     = 1 !< The reduction method is min value
+  INTEGER, PARAMETER :: time_max     = 2 !< The reduction method is max value
+  INTEGER, PARAMETER :: time_sum     = 3 !< The reduction method is sum of values
+  INTEGER, PARAMETER :: time_average= 4 !< The reduction method is average of values
+  INTEGER, PARAMETER :: time_rms     = 5 !< The reudction method is root mean square of values
   INTEGER, PARAMETER :: time_diurnal = 6 !< The reduction method is diurnal
-  INTEGER, PARAMETER :: time_power   = 7 !< The reduction method is power
+  INTEGER, PARAMETER :: time_power   = 7 !< The reduction method is average with exponents
   CHARACTER(len=7)   :: avg_name = 'average' !< Name of the average fields
   CHARACTER(len=8)   :: no_units = "NO UNITS"!< String indicating that the variable has no units
   INTEGER, PARAMETER :: begin_time  = 1 !< Use the begining of the time average bounds

diag_manager/fms_diag_axis_object.F90

@@ -371,8 +371,9 @@ subroutine write_axis_metadata(this, fms2io_fileobj, edges_in_file, parent_axis)
     end select
 
     !< Write its metadata
-    call register_variable_attribute(fms2io_fileobj, axis_name, "long_name", diag_axis%long_name, &
-      str_len=len_trim(diag_axis%long_name))
+    if(allocated(diag_axis%long_name)) &
+      call register_variable_attribute(fms2io_fileobj, axis_name, "long_name", diag_axis%long_name, &
+        str_len=len_trim(diag_axis%long_name))
 
     if (diag_axis%cart_name .NE. "N") &
       call register_variable_attribute(fms2io_fileobj, axis_name, "axis", diag_axis%cart_name, str_len=1)

diag_manager/fms_diag_field_object.F90

@@ -81,6 +81,7 @@ module fms_diag_field_object_mod
                                                                            !! the corresponding index in
                                                                            !! buffer_ids(:) is allocated.
      logical, allocatable                             :: mask(:,:,:,:)     !< Mask passed in send_data
+     logical                                          :: halo_present = .false. !< set if any halos are used
   contains
 !     procedure :: send_data => fms_send_data  !!TODO
 ! Get ID functions
@@ -168,6 +169,10 @@ module fms_diag_field_object_mod
      procedure :: get_file_ids
      procedure :: set_mask
      procedure :: allocate_mask
+     procedure :: set_halo_present
+     procedure :: is_halo_present
+     procedure :: find_missing_value
+     procedure :: has_mask_allocated
 end type fmsDiagField_type
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! variables !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 type(fmsDiagField_type) :: null_ob
@@ -1652,10 +1657,6 @@ subroutine allocate_mask(this, mask_in, omp_axis)
   class(fmsDiagAxisContainer_type), intent(in), optional :: omp_axis(:) !< true if calling from omp region
   integer :: axis_num, length(4)
   integer, pointer :: id_num
-  if(allocated(this%mask)) then
-    call mpp_error(NOTE,"set_mask:: mask already allocated for field"//this%longname)
-    deallocate(this%mask)
-  endif
   ! if not omp just allocate to whatever is given
   if(.not. present(omp_axis)) then
     allocate(this%mask(size(mask_in,1), size(mask_in,2), size(mask_in,3), &
@@ -1692,5 +1693,49 @@ subroutine set_mask(this, mask_in, is, js, ks, ie, je, ke)
   endif
 end subroutine set_mask
 
+!> sets halo_present to true
+subroutine set_halo_present(this)
+  class(fmsDiagField_type), intent(inout) :: this !< field object to modify
+  this%halo_present = .true.
+end subroutine set_halo_present
+
+!> Getter for halo_present
+pure function is_halo_present(this)
+  class(fmsDiagField_type), intent(in) :: this !< field object to get from
+  logical :: is_halo_present
+  is_halo_present = this%halo_present
+end function is_halo_present
+
+!> Helper routine to find and set the netcdf missing value for a field
+!! Always returns r8 due to reduction routine args
+!! casts up to r8 from given missing val or default if needed
+function find_missing_value(this, missing_val) &
+  result(res)
+  class(fmsDiagField_type), intent(in) :: this !< field object to get missing value for
+  class(*), allocatable, intent(out) :: missing_val !< outputted netcdf missing value (oriignal type)
+  real(r8_kind) :: res !< returned r8 copy of missing_val
+
+  if(this%has_missing_value()) then
+    missing_val = this%get_missing_value(this%get_vartype())
+  else
+    missing_val = get_default_missing_value(this%get_vartype())
+  endif
+
+  select type(missing_val)
+    type is (real(r8_kind))
+      res = missing_val
+    type is (real(r4_kind))
+      res = real(missing_val, r8_kind)
+  end select
+end function find_missing_value
+
+!> @returns allocation status of logical mask array
+!! this just indicates whether the mask array itself has been alloc'd
+!! this is different from @ref has_mask_variant, which is set earlier for whether a mask is being used at all
+pure logical function has_mask_allocated(this)
+  class(fmsDiagField_type),intent(in) :: this !< field object to check mask allocation for
+  has_mask_allocated = allocated(this%mask)
+end function has_mask_allocated
+
 #endif
 end module fms_diag_field_object_mod

diag_manager/fms_diag_file_object.F90

@@ -145,6 +145,8 @@ module fms_diag_file_object_mod
  procedure, public :: has_file_varlist
  procedure, public :: has_file_global_meta
  procedure, public :: dump_file_obj
+ procedure, public :: get_buffer_ids
+ procedure, public :: get_number_of_buffers
 end type fmsDiagFile_type
 
 type, extends (fmsDiagFile_type) :: subRegionalFile_type
@@ -1475,5 +1477,21 @@ subroutine close_diag_file(this)
   endif
 end subroutine close_diag_file
 
+!> \brief Gets the buffer_id list from the file object
+pure function get_buffer_ids (this)
+  class(fmsDiagFile_type), intent(in) :: this !< The file object
+  integer, allocatable :: get_buffer_ids(:) !< returned buffer ids for this file
+
+  allocate(get_buffer_ids(this%number_of_buffers))
+  get_buffer_ids = this%buffer_ids
+end function get_buffer_ids
+
+!> Gets the stored number of buffers from the file object
+pure function get_number_of_buffers(this)
+  class(fmsDiagFile_type), intent(in) :: this !< file object
+  integer :: get_number_of_buffers !< returned number of buffers
+  get_number_of_buffers = this%number_of_buffers
+end function get_number_of_buffers
+
 #endif
 end module fms_diag_file_object_mod

diag_manager/fms_diag_object.F90

@@ -238,6 +238,7 @@ integer function fms_register_diag_field_obj &
     bufferptr => this%FMS_diag_output_buffers(fieldptr%buffer_ids(i))
     call bufferptr%set_field_id(this%registered_variables)
     call bufferptr%set_yaml_id(fieldptr%buffer_ids(i))
+    call bufferptr%init_buffer_time(init_time)
   enddo
 
 !> Allocate and initialize member buffer_allocated of this field
@@ -538,6 +539,10 @@ logical function fms_diag_accept_data (this, diag_field_id, field_data, mask, rm
   !< Set the field_weight. If "weight" is not present it will be set to 1.0_r8_kind
   field_weight = set_weight(weight)
 
+  !< Set the variable type based off passed in field data
+  if(.not. this%FMS_diag_fields(diag_field_id)%has_vartype()) &
+    call this%FMS_diag_fields(diag_field_id)%set_type(field_data(1,1,1,1))
+
   !< Check that the indices are present in the correct combination
   error_string = check_indices_order(is_in, ie_in, js_in, je_in)
   if (trim(error_string) .ne. "") call mpp_error(FATAL, trim(error_string)//". "//trim(field_info))
@@ -550,6 +555,8 @@ logical function fms_diag_accept_data (this, diag_field_id, field_data, mask, rm
   if ((present(is_in) .and. present(ie_in)) .or. (present(js_in) .and. present(je_in))) &
     has_halos = .true.
 
+  if(has_halos) call this%FMS_diag_fields(diag_field_id)%set_halo_present()
+
   !< If the field has `mask_variant=.true.`, check that mask OR rmask are present
   if (this%FMS_diag_fields(diag_field_id)%is_mask_variant()) then
     if (.not. allocated(mask) .and. .not. allocated(rmask)) call mpp_error(FATAL, &
@@ -602,7 +609,8 @@ logical function fms_diag_accept_data (this, diag_field_id, field_data, mask, rm
     if (.not. this%FMS_diag_fields(diag_field_id)%is_data_buffer_allocated()) then
       data_buffer_is_allocated = &
         this%FMS_diag_fields(diag_field_id)%allocate_data_buffer(field_data, this%diag_axis)
-      call this%FMS_diag_fields(diag_field_id)%allocate_mask(oor_mask, this%diag_axis)
+      if(.not. this%FMS_diag_fields(diag_field_id)%has_mask_allocated()) &
+        call this%FMS_diag_fields(diag_field_id)%allocate_mask(oor_mask, this%diag_axis)
     endif
     call this%FMS_diag_fields(diag_field_id)%set_data_buffer_is_allocated(.TRUE.)
     call this%FMS_diag_fields(diag_field_id)%set_math_needs_to_be_done(.TRUE.)
@@ -621,7 +629,8 @@ logical function fms_diag_accept_data (this, diag_field_id, field_data, mask, rm
       bounds, using_blocking, Time=Time)
     if (trim(error_string) .ne. "") call mpp_error(FATAL, trim(error_string)//". "//trim(field_info))
     call this%FMS_diag_fields(diag_field_id)%set_math_needs_to_be_done(.FALSE.)
-    call this%FMS_diag_fields(diag_field_id)%allocate_mask(oor_mask)
+    if(.not. this%FMS_diag_fields(diag_field_id)%has_mask_allocated()) &
+      call this%FMS_diag_fields(diag_field_id)%allocate_mask(oor_mask)
     call this%FMS_diag_fields(diag_field_id)%set_mask(oor_mask)
     return
   end if main_if
@@ -654,7 +663,7 @@ subroutine fms_diag_send_complete(this, time_step)
   character(len=128) :: error_string
   type(fmsDiagIbounds_type) :: bounds
   integer, dimension(:), allocatable :: file_ids !< Array of file IDs for a field
-  logical, parameter :: DEBUG_SC = .true. !< turn on output for debugging
+  logical, parameter :: DEBUG_SC = .false. !< turn on output for debugging
 
   !< Update the current model time by adding the time_step
   this%current_model_time = this%current_model_time + time_step
@@ -699,18 +708,30 @@ end subroutine fms_diag_send_complete
 
 !> @brief Loops through all the files, open the file, writes out axis and
 !! variable metadata and data when necessary.
+!! TODO: passing in the saved mask from the field obj to diag_reduction_done_wrapper
+!! for performance
 subroutine fms_diag_do_io(this, is_end_of_run)
   class(fmsDiagObject_type), target, intent(inout)  :: this          !< The diag object
   logical,                 optional, intent(in)     :: is_end_of_run !< If .true. this is the end of the run,
                                                                      !! so force write
 #ifdef use_yaml
   integer :: i !< For do loops
   class(fmsDiagFileContainer_type), pointer :: diag_file !< Pointer to this%FMS_diag_files(i) (for convenience)
+  class(fmsDiagOutputBuffer_type), pointer  :: diag_buff !< pointer to output buffers iterated in buff_loop
+  class(fmsDiagField_type), pointer         :: diag_field !< pointer to output buffers iterated in buff_loop
+  class(DiagYamlFilesVar_type), pointer     :: field_yaml !< Pointer to a field from yaml fields
   TYPE (time_type),                 pointer :: model_time!< The current model time
-
+  integer, allocatable                      :: buff_ids(:) !< ids for output buffers to loop through
+  integer                                   :: ibuff !< buffer index
   logical :: file_is_opened_this_time_step !< True if the file was opened in this time_step
                                            !! If true the metadata will need to be written
-  logical :: force_write
+  logical :: force_write !< force the last write if at end of run
+  logical :: is_writing !< true if we are writing the actual field data (metadata is always written)
+  logical :: has_mask !< whether we have a mask
+  logical, parameter :: DEBUG_REDUCT = .false. !< enables debugging output
+  class(*), allocatable :: missing_val !< netcdf missing value for a given field
+  real(r8_kind) :: mval !< r8 copy of missing value
+  character(len=128) :: error_string !< outputted error string from reducti
 
   force_write = .false.
   if (present (is_end_of_run)) force_write = .true.
@@ -732,7 +753,38 @@ subroutine fms_diag_do_io(this, is_end_of_run)
       call diag_file%write_axis_data(this%diag_axis)
     endif
 
-    if (diag_file%is_time_to_write(model_time)) then
+    is_writing = diag_file%is_time_to_write(model_time)
+
+    ! finish reduction method if its time to write
+    buff_reduct: if (is_writing) then
+      allocate(buff_ids(diag_file%FMS_diag_file%get_number_of_buffers()))
+      buff_ids = diag_file%FMS_diag_file%get_buffer_ids()
+      ! loop through the buffers and finish reduction if needed
+      buff_loop: do ibuff=1, SIZE(buff_ids)
+        diag_buff => this%FMS_diag_output_buffers(buff_ids(ibuff))
+        field_yaml => diag_yaml%get_diag_field_from_id(diag_buff%get_yaml_id())
+        diag_field => this%FMS_diag_fields(diag_buff%get_field_id())
+        ! sets missing value
+        mval = diag_field%find_missing_value(missing_val)
+        ! time_average and greater values all involve averaging so need to be "finished" before written
+        if( field_yaml%has_var_reduction()) then
+          if( field_yaml%get_var_reduction() .ge. time_average) then
+            if(DEBUG_REDUCT)call mpp_error(NOTE, "fms_diag_do_io:: finishing reduction for "//diag_field%get_longname())
+            has_mask = diag_field%has_mask_variant()
+            if(has_mask) has_mask = diag_field%get_mask_variant()
+            error_string = diag_buff%diag_reduction_done_wrapper( &
+                                    field_yaml%get_var_reduction(), &
+                                    mval, has_mask)
+          endif
+        endif
+        !endif
+        nullify(diag_buff)
+        nullify(field_yaml)
+      enddo buff_loop
+      deallocate(buff_ids)
+    endif buff_reduct
+
+    if (is_writing) then
       call diag_file%increase_unlim_dimension_level()
       call diag_file%write_time_data()
       call diag_file%write_field_data(this%FMS_diag_fields, this%FMS_diag_output_buffers)
@@ -795,6 +847,8 @@ function fms_diag_do_reduction(this, field_data, diag_field_id, oor_mask, weight
   real(kind=r8_kind)        :: missing_value      !< Missing_value for data points that are masked
                                                   !! This will obtained as r8 and converted to the right type as
                                                   !! needed. This is to avoid yet another select type ...
+  logical                   :: new_time           !< .True. if this is a new time (i.e data has not be been
+                                                  !! sent for this time)
 
   !TODO mostly everything
   field_ptr => this%FMS_diag_fields(diag_field_id)
@@ -908,11 +962,17 @@ function fms_diag_do_reduction(this, field_data, diag_field_id, oor_mask, weight
       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)
+        bounds_in, bounds_out, missing_value, .true.)
       if (trim(error_msg) .ne. "") then
         return
       endif
     case (time_average)
+      new_time = buffer_ptr%update_buffer_time(time)
+      error_msg = buffer_ptr%do_time_sum_wrapper(field_data, oor_mask, field_ptr%get_mask_variant(), &
+        bounds_in, bounds_out, missing_value, new_time)
+      if (trim(error_msg) .ne. "") then
+        return
+      endif
     case (time_power)
     case (time_rms)
     case (time_diurnal)

diag_manager/fms_diag_output_buffer.F90

@@ -27,14 +27,14 @@ module fms_diag_output_buffer_mod
 #ifdef use_yaml
 use platform_mod
 use iso_c_binding
-use time_manager_mod, only: time_type, operator(==)
-use mpp_mod, only: mpp_error, FATAL
+use time_manager_mod, only: time_type, operator(==), operator(>)
+use mpp_mod, only: mpp_error, FATAL, NOTE
 use diag_data_mod, only: DIAG_NULL, DIAG_NOT_REGISTERED, i4, i8, r4, r8, get_base_time, MIN_VALUE, MAX_VALUE, EMPTY, &
                          time_min, time_max
 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, do_time_sum_update
+use fms_diag_reduction_methods_mod, only: do_time_none, do_time_min, do_time_max, do_time_sum_update, time_update_done
 use fms_diag_time_utils_mod, only: diag_time_inc
 
 implicit none
@@ -54,6 +54,7 @@ module fms_diag_output_buffer_mod
   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
   logical               :: done_with_math     !< .True. if done doing the math
+  type(time_type)       :: time               !< The last time the data was received
 
   contains
   procedure :: add_axis_ids
@@ -62,6 +63,8 @@ module fms_diag_output_buffer_mod
   procedure :: get_field_id
   procedure :: set_yaml_id
   procedure :: get_yaml_id
+  procedure :: init_buffer_time
+  procedure :: update_buffer_time
   procedure :: is_done_with_math
   procedure :: set_done_with_math
   procedure :: write_buffer
@@ -77,7 +80,8 @@ module fms_diag_output_buffer_mod
   procedure :: do_time_min_wrapper
   procedure :: do_time_max_wrapper
   procedure :: do_time_sum_wrapper
-
+  procedure :: diag_reduction_done_wrapper
+  procedure :: get_buffer_dims
 end type fmsDiagOutputBuffer_type
 
 ! public types
@@ -323,6 +327,35 @@ subroutine set_yaml_id(this, yaml_id)
   this%yaml_id = yaml_id
 end subroutine set_yaml_id
 
+!> @brief inits the buffer time for the buffer
+subroutine init_buffer_time(this, time)
+  class(fmsDiagOutputBuffer_type), intent(inout) :: this        !< Buffer object
+  type(time_type), optional,       intent(in)    :: time        !< time to add to the buffer
+
+  if (present(time)) then
+    this%time = time
+  else
+    this%time = get_base_time()
+  endif
+end subroutine init_buffer_time
+
+!> @brief Update the buffer time if it is a new time
+!! @return .true. if the buffer was updated
+function update_buffer_time(this, time) &
+  result(res)
+  class(fmsDiagOutputBuffer_type), intent(inout) :: this        !< Buffer object
+  type(time_type),                 intent(in)    :: time        !< time to add to the buffer
+
+  logical :: res
+
+  if (time > this%time) then
+    this%time = time
+    res = .true.
+  else
+    res = .false.
+  endif
+end function
+
 !> @brief Determine if finished with math
 !! @return this%done_with_math
 function is_done_with_math(this) &
@@ -554,7 +587,8 @@ 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) &
+function do_time_sum_wrapper(this, field_data, mask, is_masked, bounds_in, bounds_out, missing_value, &
+                             increase_counter) &
   result(err_msg)
   class(fmsDiagOutputBuffer_type), intent(inout) :: this                !< buffer object to write
   class(*),                        intent(in)    :: field_data(:,:,:,:) !< Buffer data for current time
@@ -563,6 +597,8 @@ function do_time_sum_wrapper(this, field_data, mask, is_masked, bounds_in, bound
   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
+  logical,                         intent(in)    :: increase_counter    !< .True. if data has not been received for
+                                                                        !! time, so the counter needs to be increased
   character(len=50) :: err_msg
 
   !TODO This will be expanded for integers
@@ -572,21 +608,55 @@ function do_time_sum_wrapper(this, field_data, mask, is_masked, bounds_in, bound
       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)
+                                bounds_in, bounds_out, missing_value, increase_counter)
       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))
+          real(missing_value, kind=r4_kind), increase_counter)
       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
+
+!> Finishes calculations for any reductions that use an average (avg, rms, pow)
+!! TODO add mask and any other needed args for adjustment, and pass in the adjusted mask
+!! to time_update_done
+function diag_reduction_done_wrapper(this, reduction_method, missing_value, has_mask) & !! , has_halo, mask) &
+  result(err_msg)
+  class(fmsDiagOutputBuffer_type), intent(inout) :: this !< Updated buffer object
+  integer, intent(in)                            :: reduction_method !< enumerated reduction type from diag_data
+  real(kind=r8_kind), intent(in)                 :: missing_value !< missing_value for masked data points
+  logical, intent(in)                            :: has_mask !< indicates if there was a mask used during buffer updates
+  character(len=51)                              :: err_msg !< error message to return, blank if sucessful
+
+  if(.not. allocated(this%buffer)) return
+
+  if(this%weight_sum .eq. 0.0_r8_kind) return
+
+  err_msg = ""
+  select type(buff => this%buffer)
+    type is (real(r8_kind))
+      call time_update_done(buff, this%weight_sum, reduction_method, missing_value, has_mask)
+    type is (real(r4_kind))
+      call time_update_done(buff, this%weight_sum, reduction_method, real(missing_value, r4_kind), has_mask)
+  end select
+  this%weight_sum = 0.0_r8_kind
+
+end function
+
+!> this leaves out the diurnal index cause its only used for tmp mask allocation
+pure function get_buffer_dims(this)
+  class(fmsDiagOutputBuffer_type), intent(in) :: this
+  integer :: get_buffer_dims(4)
+  get_buffer_dims = this%buffer_dims(1:4)
+end function
+
 #endif
 end module fms_diag_output_buffer_mod

diag_manager/fms_diag_reduction_methods.F90

@@ -30,12 +30,13 @@
 module fms_diag_reduction_methods_mod
   use platform_mod, only: r8_kind, r4_kind
   use fms_diag_bbox_mod, only: fmsDiagIbounds_type
+  use fms_string_utils_mod, only: string
   use mpp_mod
   implicit none
   private
 
   public :: check_indices_order, init_mask, set_weight
-  public :: do_time_none, do_time_min, do_time_max, do_time_sum_update
+  public :: do_time_none, do_time_min, do_time_max, do_time_sum_update, time_update_done
 
   !> @brief Does the time_none reduction method. See include/fms_diag_reduction_methods.inc
   !TODO This needs to be extended to integers
@@ -62,6 +63,13 @@ module fms_diag_reduction_methods_mod
     module procedure do_time_sum_update_r4, do_time_sum_update_r8
   end interface
 
+  !> @brief Finishes a reduction that involves an average
+  !! (ie. time_avg, rms, pow)
+  !! This takes the average at the end of the time step
+  interface time_update_done
+    module procedure sum_update_done_r4, sum_update_done_r8
+  end interface
+
   contains
 
   !> @brief Checks improper combinations of is, ie, js, and je.

diag_manager/include/fms_diag_reduction_methods.inc

@@ -19,7 +19,7 @@
 
 ! for any debug prints
 #ifndef DEBUG_REDUCT
-#define DEBUG_REDUCT .true.
+#define DEBUG_REDUCT .false.
 #endif
 
 !> @brief Do the time_none reduction method (i.e copy the correct portion of the input data)
@@ -215,7 +215,7 @@ end subroutine DO_TIME_MAX_
 !!
 !! Where l are the indices passed in through the bounds_in/out
 subroutine DO_TIME_SUM_UPDATE_(data_out, weight_sum, data_in, mask, is_masked, bounds_in, bounds_out, &
-                               missing_value, weight, pow)
+                               missing_value, increase_counter, weight, pow)
   real(FMS_TRM_KIND_),       intent(inout) :: data_out(:,:,:,:,:) !< output data
   real(r8_kind),       intent(inout) :: weight_sum          !< Sum of weights from the output buffer object
   real(FMS_TRM_KIND_),       intent(in)    :: data_in(:,:,:,:)    !< data to update the buffer with
@@ -226,6 +226,8 @@ subroutine DO_TIME_SUM_UPDATE_(data_out, weight_sum, data_in, mask, is_masked, b
   type(fmsDiagIbounds_type), intent(in)    :: bounds_out          !< indices indicating the correct portion
                                                                   !! of the output buffer
   real(FMS_TRM_KIND_),       intent(in)    :: missing_value       !< Missing_value for data points that are masked
+  logical,                   intent(in)    :: increase_counter    !< .True. if data has not been received for
+                                                                  !! time, so the counter needs to be increased
   real(r8_kind),optional, intent(in) :: weight              !< Weight applied to data_in before added to data_out
                                                                   !! used for weighted averages, default 1.0
   real(FMS_TRM_KIND_),optional, intent(in) :: pow                 !< Used for pow reduction, adds field^pow to buffer
@@ -251,7 +253,7 @@ subroutine DO_TIME_SUM_UPDATE_(data_out, weight_sum, data_in, mask, is_masked, b
   endif
 
   ! update with given weight for average before write
-  weight_sum = weight_sum + weight_loc
+  if (increase_counter) weight_sum = weight_sum + weight_loc
 
   is_out = bounds_out%get_imin()
   ie_out = bounds_out%get_imax()
@@ -289,11 +291,38 @@ subroutine DO_TIME_SUM_UPDATE_(data_out, weight_sum, data_in, mask, is_masked, b
     do k = 0, ke_out - ks_out
       do j = 0, je_out - js_out
         do i = 0, ie_out - is_out
-          data_out(is_out + i, js_out + j, ks_out + k,  :, 1) =   &
+          data_out(is_out + i, js_out + j, ks_out + k,  :, 1) =  &
              data_out(is_out + i, js_out + j, ks_out + k, :, 1)  &
             + (data_in(is_in +i, js_in + j, ks_in + k, :) * weight_loc) ** pow_loc
         enddo
       enddo
     enddo
   endif
 end subroutine DO_TIME_SUM_UPDATE_
+
+!> To be called with diag_send_complete, finishes reductions
+!! Just divides the buffer by the counter array(which is just the sum of the weights used in the buffer's reduction)
+!! TODO: change has_mask to an actual logical mask so we don't have to check for missing values
+subroutine SUM_UPDATE_DONE_(out_buffer_data, weight_sum, reduction_method, missing_val, has_mask)
+  real(FMS_TRM_KIND_), intent(inout) :: out_buffer_data(:,:,:,:,:) !< data buffer previously updated with
+                                                                   !! do_time_sum_update
+  real(r8_kind), intent(in)          :: weight_sum !< sum of weights for averaging, provided via argument to send data
+  integer, intent(in)                :: reduction_method !< which reduction method to use, should be time_avg
+  real(FMS_TRM_KIND_), intent(in)    :: missing_val !< missing value for masked elements
+  logical, intent(in)                :: has_mask !< indicates if mask is used so missing values can be skipped
+  !! TODO replace conditional in the `where` with passed in and ajusted mask from the original call
+  !logical, optional, intent(in)      :: mask(:,:,:,:) !< logical mask from accept data call, if using one.
+  !logical                            :: has_mask !< whether or not mask is present
+
+  if ( has_mask ) then
+    where(out_buffer_data(:,:,:,:,1) .ne. missing_val)
+      out_buffer_data(:,:,:,:,1) = out_buffer_data(:,:,:,:,1) &
+                                   / weight_sum
+    endwhere
+  else  !not mask variant
+    out_buffer_data(:,:,:,:,1) = out_buffer_data(:,:,:,:,1) &
+                               / weight_sum
+  endif
+
+end subroutine
+

diag_manager/include/fms_diag_reduction_methods_r4.fh

@@ -38,6 +38,9 @@
 #undef DO_TIME_SUM_UPDATE_
 #define DO_TIME_SUM_UPDATE_ do_time_sum_update_r4
 
+#undef SUM_UPDATE_DONE_
+#define SUM_UPDATE_DONE_ sum_update_done_r4
+
 #include "fms_diag_reduction_methods.inc"
 
 !> @}

diag_manager/include/fms_diag_reduction_methods_r8.fh

@@ -38,6 +38,9 @@
 #undef DO_TIME_SUM_UPDATE_
 #define DO_TIME_SUM_UPDATE_ do_time_sum_update_r8
 
+#undef SUM_UPDATE_DONE_
+#define SUM_UPDATE_DONE_ sum_update_done_r8
+
 #include "fms_diag_reduction_methods.inc"
 
 !> @}

test_fms/diag_manager/Makefile.am

@@ -31,7 +31,7 @@ LDADD = $(top_builddir)/libFMS/libFMS.la
 check_PROGRAMS = test_diag_manager test_diag_manager_time \
  test_diag_dlinked_list test_diag_yaml test_diag_ocean test_modern_diag test_diag_buffer \
  test_flexible_time test_diag_update_buffer test_reduction_methods check_time_none \
- check_time_min check_time_max check_time_sum
+ check_time_min check_time_max check_time_sum check_time_avg
 
 # This is the source code for the test.
 test_diag_manager_SOURCES = test_diag_manager.F90
@@ -48,19 +48,21 @@ check_time_none_SOURCES = testing_utils.F90 check_time_none.F90
 check_time_min_SOURCES = testing_utils.F90 check_time_min.F90
 check_time_max_SOURCES = testing_utils.F90 check_time_max.F90
 check_time_sum_SOURCES = testing_utils.F90 check_time_sum.F90
+check_time_avg_SOURCES = testing_utils.F90 check_time_avg.F90
 
 TEST_EXTENSIONS = .sh
 SH_LOG_DRIVER = env AM_TAP_AWK='$(AWK)' $(SHELL) \
                   $(abs_top_srcdir)/test_fms/tap-driver.sh
 
 # Run the test.
-TESTS = test_diag_manager2.sh test_time_none.sh test_time_min.sh test_time_max.sh test_time_sum.sh
+TESTS = test_diag_manager2.sh test_time_none.sh test_time_min.sh test_time_max.sh test_time_sum.sh \
+        test_time_avg.sh
 
 testing_utils.mod: testing_utils.$(OBJEXT)
 
 # Copy over other needed files to the srcdir
 EXTRA_DIST = test_diag_manager2.sh check_crashes.sh test_time_none.sh test_time_min.sh test_time_max.sh \
-             test_time_sum.sh
+             test_time_sum.sh test_time_avg.sh
 
 if USING_YAML
 skipflag=""

test_fms/diag_manager/check_time_avg.F90

@@ -0,0 +1,270 @@
+!***********************************************************************
+!*                   GNU Lesser General Public License
+!*
+!* This file is part of the GFDL Flexible Modeling System (FMS).
+!*
+!* FMS is free software: you can redistribute it and/or modify it under
+!* the terms of the GNU Lesser General Public License as published by
+!* the Free Software Foundation, either version 3 of the License, or (at
+!* your option) any later version.
+!*
+!* FMS is distributed in the hope that it will be useful, but WITHOUT
+!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+!* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+!* for more details.
+!*
+!* You should have received a copy of the GNU Lesser General Public
+!* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
+!***********************************************************************
+
+!> @brief  Checks the output file after running test_reduction_methods using the "time_avg" reduction method
+program check_time_avg
+  use fms_mod,           only: fms_init, fms_end, string
+  use fms2_io_mod,       only: FmsNetcdfFile_t, read_data, close_file, open_file
+  use mpp_mod,           only: mpp_npes, mpp_error, FATAL, mpp_pe, input_nml_file
+  use platform_mod,      only: r4_kind, r8_kind
+  use testing_utils,     only: allocate_buffer, test_normal, test_openmp, test_halos, no_mask, logical_mask, real_mask
+
+  implicit none
+
+  type(FmsNetcdfFile_t)              :: fileobj            !< FMS2 fileobj
+  type(FmsNetcdfFile_t)              :: fileobj1           !< FMS2 fileobj for subregional file 1
+  type(FmsNetcdfFile_t)              :: fileobj2           !< FMS2 fileobj for subregional file 2
+  real(kind=r4_kind), allocatable    :: cdata_out(:,:,:,:) !< Data in the compute domain
+  integer                            :: nx                 !< Number of points in the x direction
+  integer                            :: ny                 !< Number of points in the y direction
+  integer                            :: nz                 !< Number of points in the z direction
+  integer                            :: nw                 !< Number of points in the 4th dimension
+  integer                            :: ti                  !< For looping through time levels
+  integer                            :: io_status          !< Io status after reading the namelist
+  logical                            :: use_mask           !< .true. if using masks
+  integer, parameter :: file_freq = 6 !< file frequency as set in diag_table.yaml
+
+  integer :: test_case = test_normal !< Indicates which test case to run
+  integer :: mask_case = no_mask     !< Indicates which masking option to run
+  integer, parameter :: kindl = KIND(0.0) !< compile-time default kind size
+
+  namelist / test_reduction_methods_nml / test_case, mask_case
+
+  call fms_init()
+
+  read (input_nml_file, test_reduction_methods_nml, iostat=io_status)
+
+  select case(mask_case)
+  case (no_mask)
+    use_mask = .false.
+  case (logical_mask, real_mask)
+    use_mask = .true.
+  end select
+  nx = 96
+  ny = 96
+  nz = 5
+  nw = 2
+
+  if (.not. open_file(fileobj, "test_avg.nc", "read")) &
+    call mpp_error(FATAL, "unable to open test_avg.nc")
+
+  if (.not. open_file(fileobj1, "test_avg_regional.nc.0004", "read")) &
+    call mpp_error(FATAL, "unable to open test_avg_regional.nc.0004")
+
+  if (.not. open_file(fileobj2, "test_avg_regional.nc.0005", "read")) &
+    call mpp_error(FATAL, "unable to open test_avg_regional.nc.0005")
+
+  cdata_out = allocate_buffer(1, nx, 1, ny, nz, nw)
+
+  do ti = 1, 8
+    cdata_out = -999_r4_kind
+    print *, "Checking answers for var0_avg - time_level:", string(ti)
+    call read_data(fileobj, "var0_avg", cdata_out(1,1,1,1), unlim_dim_level=ti)
+    call check_data_0d(cdata_out(1,1,1,1), ti)
+
+    cdata_out = -999_r4_kind
+    print *, "Checking answers for var1_avg - time_level:", string(ti)
+    call read_data(fileobj, "var1_avg", cdata_out(:,1,1,1), unlim_dim_level=ti)
+    call check_data_1d(cdata_out(:,1,1,1), ti)
+
+    cdata_out = -999_r4_kind
+    print *, "Checking answers for var2_avg - time_level:", string(ti)
+    call read_data(fileobj, "var2_avg", cdata_out(:,:,1,1), unlim_dim_level=ti)
+    call check_data_2d(cdata_out(:,:,1,1), ti)
+
+    cdata_out = -999_r4_kind
+    print *, "Checking answers for var3_avg - time_level:", string(ti)
+    call read_data(fileobj, "var3_avg", cdata_out(:,:,:,1), unlim_dim_level=ti)
+    call check_data_3d(cdata_out(:,:,:,1), ti, .false.)
+
+    cdata_out = -999_r4_kind
+    print *, "Checking answers for var4_avg - time_level:", string(ti)
+    call read_data(fileobj, "var4_avg", cdata_out(:,:,:,:), unlim_dim_level=ti)
+    call check_data_3d(cdata_out(:,:,:,1), ti, .false.)
+    call check_data_3d(cdata_out(:,:,:,2), ti, .false.)
+
+    cdata_out = -999_r4_kind
+    print *, "Checking answers for var3_Z - time_level:", string(ti)
+    call read_data(fileobj, "var3_Z", cdata_out(:,:,1:2,1), unlim_dim_level=ti)
+    call check_data_3d(cdata_out(:,:,1:2,1), ti, .true., nz_offset=1)
+
+    cdata_out = -999_r4_kind
+    print *, "Checking answers for var3_avg in the first regional file- time_level:", string(ti)
+    call read_data(fileobj1, "var3_avg", cdata_out(1:4,1:3,1:2,1), unlim_dim_level=ti)
+    call check_data_3d(cdata_out(1:4,1:3,1:2,1), ti, .true., nx_offset=77, ny_offset=77, nz_offset=1)
+
+    cdata_out = -999_r4_kind
+    print *, "Checking answers for var3_avg in the second regional file- time_level:", string(ti)
+    call read_data(fileobj2, "var3_avg", cdata_out(1:4,1:1,1:2,1), unlim_dim_level=ti)
+    call check_data_3d(cdata_out(1:4,1:1,1:2,1), ti, .true., nx_offset=77, ny_offset=80, nz_offset=1)
+  enddo
+
+  call fms_end()
+
+contains
+
+  ! sent data set to:
+  !          buffer(ii-is+1+nhalo, j-js+1+nhalo, k, l) = real(ii, kind=r8_kind)* 1000_r8_kind + &
+  !                                                      real(j, kind=r8_kind)* 10_r8_kind + &
+  !                                                      real(k, kind=r8_kind)
+  !                                                      + time_index/100
+  !> @brief Check that the 0d data read in is correct
+  subroutine check_data_0d(buffer, time_level)
+    real(kind=r4_kind), intent(inout) :: buffer        !< Buffer read from the table
+    integer,            intent(in)    :: time_level    !< Time level read in
+
+    real(kind=r4_kind)                :: buffer_exp    !< Expected result
+    integer :: i, step_avg = 0 !< avg of time step increments to use in generating reference data
+
+    ! avgs integers for decimal part of field input
+    ! ie. level 1 = 1+2+..+6
+    !           2 = 7+8+..+12
+    step_avg = 0
+    do i=(time_level-1)*file_freq+1, time_level*file_freq
+      step_avg = step_avg + i
+    enddo
+
+   ! 0d answer is:
+   !    (1011 * frequency avg'd over )
+   !  + ( 1/100 * avg of time step increments )
+    buffer_exp = real((1000.0_r8_kind+10.0_r8_kind+1.0_r8_kind) * file_freq + &
+                      real(step_avg,r8_kind)/100.0_r8_kind, kind=r4_kind)
+    buffer_exp = buffer_exp / file_freq
+
+    if (abs(buffer - buffer_exp) > 0.0) print *, "answer not exact for 0d, time:", time_level, &
+                                                 " diff:", abs(buffer-buffer_exp)
+
+    if (abs(buffer - buffer_exp) > 1.0e-4) then
+      print *, "time_level", time_level, "expected", buffer_exp, "read", buffer
+      call mpp_error(FATAL, "Check_time_avg::check_data_0d:: Data is not correct")
+    endif
+  end subroutine check_data_0d
+
+  !> @brief Check that the 1d data read in is correct
+  subroutine check_data_1d(buffer, time_level)
+    real(kind=r4_kind), intent(in)    :: buffer(:)     !< Buffer read from the table
+    integer,            intent(in)    :: time_level    !< Time level read in
+    real(kind=r4_kind)                :: buffer_exp    !< Expected result
+    integer :: step_sum !< avg of time step increments to use in generating reference data
+    integer :: ii, i, j, k, l !< For looping
+    integer :: n
+
+    step_sum = 0
+    do i=(time_level-1)*file_freq+1, time_level*file_freq
+      step_sum = step_sum + i
+    enddo
+
+    ! 1d answer is
+    !  (((i * 1000 + 11) * frequency) + (sum of time steps)) / frequency
+    !  or
+    !  => (i * 1000 + 11) + (sum of time_steps/frequency/100)
+    do ii = 1, size(buffer, 1)
+      buffer_exp =  real( &
+                        (real(ii, kind=r8_kind)*1000.0_r8_kind +11.0_r8_kind) + &
+                        (real(step_sum, kind=r8_kind)/file_freq/100.0_r8_kind) &
+                    , kind=r4_kind)
+
+      if (use_mask .and. ii .eq. 1) buffer_exp = -666_r4_kind
+      if (abs(buffer(ii) - buffer_exp) > 0.0) then
+        print *, "i:", ii, "read in:", buffer(ii), "expected:", buffer_exp, "time level:", time_level
+        print *, "diff:", abs(buffer(ii) - buffer_exp)
+        call mpp_error(FATAL, "Check_time_avg::check_data_1d:: Data is not correct")
+      endif
+    enddo
+  end subroutine check_data_1d
+
+  !> @brief Check that the 2d data read in is correct
+  subroutine check_data_2d(buffer, time_level)
+    real(kind=r4_kind), intent(in)    :: buffer(:,:)   !< Buffer read from the table
+    integer,            intent(in)    :: time_level    !< Time level read in
+    real(kind=r4_kind)                :: buffer_exp    !< Expected result
+
+    integer :: ii,i, j, k, l !< For looping
+    integer :: step_avg !< avg of time step increments to use in generating reference data
+
+    step_avg = 0
+    do i=(time_level-1)*file_freq+1, time_level*file_freq
+      step_avg = step_avg + i
+    enddo
+
+    ! 2d answer is
+    !  ((i * 1000 + j * 10 + 1) * frequency) + (avg of time steps)
+    do ii = 1, size(buffer, 1)
+      do j = 1, size(buffer, 2)
+        buffer_exp = real(real(ii, kind=r8_kind)* 1000.0_kindl+ &
+                     10.0_kindl*real(j, kind=r8_kind)+1.0_kindl + &
+                     real(step_avg, kind=r8_kind)/file_freq/100.0_r8_kind, kind=r4_kind)
+        if (use_mask .and. ii .eq. 1 .and. j .eq. 1) buffer_exp = -666_r4_kind
+        if (abs(buffer(ii, j) - buffer_exp) > 0.0) then
+          print *, "indices:", ii, j, "expected:", buffer_exp, "read in:",buffer(ii, j)
+          call mpp_error(FATAL, "Check_time_avg::check_data_2d:: Data is not correct")
+        endif
+      enddo
+    enddo
+  end subroutine check_data_2d
+
+  !> @brief Check that the 3d data read in is correct
+  subroutine check_data_3d(buffer, time_level, is_regional, nx_offset, ny_offset, nz_offset)
+    real(kind=r4_kind), intent(in)    :: buffer(:,:,:) !< Buffer read from the table
+    integer,            intent(in)    :: time_level    !< Time level read in
+    logical,            intent(in)    :: is_regional   !< .True. if the variable is subregional
+    real(kind=r4_kind)                :: buffer_exp    !< Expected result
+    integer, optional,  intent(in)    :: nx_offset     !< Offset in the x direction
+    integer, optional,  intent(in)    :: ny_offset     !< Offset in the y direction
+    integer, optional,  intent(in)    :: nz_offset     !< Offset in the z direction
+
+    integer :: ii, i, j, k, l !< For looping
+    integer :: nx_oset !< Offset in the x direction (local variable)
+    integer :: ny_oset !< Offset in the y direction (local variable)
+    integer :: nz_oset !< Offset in the z direction (local variable)
+    integer :: step_avg!< avg of time step increments to use in generating reference data
+
+    step_avg = 0
+    do i=(time_level-1)*file_freq+1, time_level*file_freq
+      step_avg = step_avg + i
+    enddo
+
+    nx_oset = 0
+    if (present(nx_offset)) nx_oset = nx_offset
+
+    ny_oset = 0
+    if (present(ny_offset)) ny_oset = ny_offset
+
+    nz_oset = 0
+    if (present(nz_offset)) nz_oset = nz_offset
+
+    ! 3d answer is
+    !  ((i * 1000 + j * 10 + k) * frequency) + (avg of time steps)
+    do ii = 1, size(buffer, 1)
+      do j = 1, size(buffer, 2)
+        do k = 1, size(buffer, 3)
+          buffer_exp = real(real(ii+nx_oset, kind=r8_kind)* 1000.0_kindl + &
+                       10.0_kindl*real(j+ny_oset, kind=r8_kind) + &
+                       1.0_kindl*real(k+nz_oset, kind=r8_kind) + &
+                       real(step_avg, kind=r8_kind)/file_freq/100.0_kindl, kind=r4_kind)
+          if (use_mask .and. ii .eq. 1 .and. j .eq. 1 .and. k .eq. 1 .and. .not. is_regional) buffer_exp = -666_r4_kind
+          if (abs(buffer(ii, j, k) - buffer_exp) > 0.0) then
+            print *, mpp_pe(),'indices:',ii, j, k, "read in:", buffer(ii, j, k), "expected:",buffer_exp
+            call mpp_error(FATAL, "Check_time_avg::check_data_3d:: Data is not correct")
+          endif
+        enddo
+      enddo
+    enddo
+  end subroutine check_data_3d
+end program

test_fms/diag_manager/test_diag_manager2.sh

@@ -713,7 +713,7 @@ diag_files:
     var_name: var5
     reduction: average
     kind: r4
-  - module: lnd_mod
+  - module: atm_mod
     var_name: var7
     reduction: average
     kind: r4

test_fms/diag_manager/test_modern_diag.F90

@@ -248,7 +248,7 @@ subroutine allocate_dummy_data(var, lat_lon_domain, cube_sphere, lnd_domain, nz)
   allocate(var%var4(is:ie, js:je, nz)) !< Variable in a 3D cube sphere domain
 
   call mpp_get_UG_compute_domain(lnd_domain, size=nland)
-  allocate(var%var5(nz)) !< Variable in the land unstructured domain
+  allocate(var%var5(nland)) !< Variable in the land unstructured domain
 
   allocate(var%var6(nz)) !< 1D variable not domain decomposed
 

test_fms/diag_manager/test_time_avg.sh

@@ -0,0 +1,180 @@
+#!/bin/sh
+
+#***********************************************************************
+#*                   GNU Lesser General Public License
+#*
+#* This file is part of the GFDL Flexible Modeling System (FMS).
+#*
+#* FMS is free software: you can redistribute it and/or modify it under
+#* the terms of the GNU Lesser General Public License as published by
+#* the Free Software Foundation, either version 3 of the License, or (at
+#* your option) any later version.
+#*
+#* FMS is distributed in the hope that it will be useful, but WITHOUT
+#* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+#* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+#* for more details.
+#*
+#* You should have received a copy of the GNU Lesser General Public
+#* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
+#***********************************************************************
+
+# Set common test settings.
+. ../test-lib.sh
+
+if [ -z "${skipflag}" ]; then
+# create and enter directory for in/output files
+output_dir
+
+cat <<_EOF > diag_table.yaml
+title: test_avg
+base_date: 2 1 1 0 0 0
+diag_files:
+- file_name: test_avg
+  time_units: hours
+  unlimdim: time
+  freq: 6 hours
+  varlist:
+  - module: ocn_mod
+    var_name: var0
+    output_name: var0_avg
+    reduction: average
+    kind: r4
+  - module: ocn_mod
+    var_name: var1
+    output_name: var1_avg
+    reduction: average
+    kind: r4
+  - module: ocn_mod
+    var_name: var2
+    output_name: var2_avg
+    reduction: average
+    kind: r4
+  - module: ocn_mod
+    var_name: var3
+    output_name: var3_avg
+    reduction: average
+    kind: r4
+  - module: ocn_mod
+    var_name: var4
+    output_name: var4_avg
+    reduction: average
+    kind: r4
+  - module: ocn_mod
+    var_name: var3
+    output_name: var3_Z
+    reduction: average
+    zbounds: 2. 3.
+    kind: r4
+- file_name: test_avg_regional
+  time_units: hours
+  unlimdim: time
+  sub_region:
+  - grid_type: latlon
+    corner1: 78. 78.
+    corner2: 78. 78.
+    corner3: 81. 81.
+    corner4: 81. 81.
+  freq: 6 hours
+  varlist:
+  - module: ocn_mod
+    var_name: var3
+    output_name: var3_avg
+    reduction: average
+    zbounds: 2. 3.
+    kind: r4
+_EOF
+
+# remove any existing files that would result in false passes during checks
+rm -f *.nc
+
+# tests with no mask, no openmp
+my_test_count=1
+printf "&diag_manager_nml \n use_modern_diag=.true. \n / \n&test_reduction_methods_nml \n test_case = 0 \n/" | cat > input.nml
+test_expect_success "Running diag_manager with "avg" reduction method (test $my_test_count)" '
+  mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "avg" reduction method (test $my_test_count)" '
+  mpirun -n 1 ../check_time_avg
+'
+
+my_test_count=`expr $my_test_count + 1`
+printf "&diag_manager_nml \n use_modern_diag=.true. \n / \n &test_reduction_methods_nml \n test_case = 0 \n mask_case = 1 \n \n/" | cat > input.nml
+test_expect_success "Running diag_manager with "avg" reduction method, logical mask (test $my_test_count)" '
+  mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "avg" reduction method, logical mask (test $my_test_count)" '
+  mpirun -n 1 ../check_time_avg
+'
+
+my_test_count=`expr $my_test_count + 1`
+printf "&diag_manager_nml \n use_modern_diag=.true. \n / \n &test_reduction_methods_nml \n test_case = 0 \n mask_case = 2 \n \n/" | cat > input.nml
+test_expect_success "Running diag_manager with "avg" reduction method, real mask (test $my_test_count)" '
+  mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "avg" reduction method, real mask (test $my_test_count)" '
+  mpirun -n 1 ../check_time_avg
+'
+
+# openmp tests
+
+export OMP_NUM_THREADS=2
+my_test_count=`expr $my_test_count + 1`
+printf "&diag_manager_nml \n use_modern_diag=.true. \n / \n&test_reduction_methods_nml \n test_case = 1 \n \n/" | cat > input.nml
+test_expect_success "Running diag_manager with "avg" reduction method with openmp (test $my_test_count)" '
+  mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "avg" reduction method with openmp (test $my_test_count)" '
+  mpirun -n 1 ../check_time_avg
+'
+
+my_test_count=`expr $my_test_count + 1`
+printf "&diag_manager_nml \n use_modern_diag=.true. \n / \n&test_reduction_methods_nml \n test_case = 1 \n mask_case = 1 \n \n/" | cat > input.nml
+test_expect_success "Running diag_manager with "avg" reduction method with openmp, logical mask (test $my_test_count)" '
+  mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "avg" reduction method with openmp, logical mask (test $my_test_count)" '
+  mpirun -n 1 ../check_time_avg
+'
+
+my_test_count=`expr $my_test_count + 1`
+printf "&diag_manager_nml \n use_modern_diag=.true. \n / \n&test_reduction_methods_nml \n test_case = 1 \n mask_case = 2 \n \n/" | cat > input.nml
+test_expect_success "Running diag_manager with "avg" reduction method with openmp, real mask (test $my_test_count)" '
+  mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "avg" reduction method with openmp, real mask (test $my_test_count)" '
+  mpirun -n 1 ../check_time_avg
+'
+
+# halo output and mask tests
+
+export OMP_NUM_THREADS=1
+
+my_test_count=`expr $my_test_count + 1`
+printf "&diag_manager_nml \n use_modern_diag=.true. \n / \n&test_reduction_methods_nml \n test_case = 2 \n \n/" | cat > input.nml
+test_expect_success "Running diag_manager with "avg" reduction method with halo output (test $my_test_count)" '
+  mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "avg" reduction method with halo output (test $my_test_count)" '
+  mpirun -n 1 ../check_time_avg
+'
+
+my_test_count=`expr $my_test_count + 1`
+printf "&diag_manager_nml \n use_modern_diag=.true. \n / \n&test_reduction_methods_nml \n test_case = 2 \n mask_case = 1 \n \n/" | cat > input.nml
+test_expect_success "Running diag_manager with "avg" reduction method with halo output with logical mask (test $my_test_count)" '
+  mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "avg" reduction method with halo output with logical mask (test $my_test_count)" '
+  mpirun -n 1 ../check_time_avg
+'
+
+my_test_count=`expr $my_test_count + 1`
+printf "&diag_manager_nml \n use_modern_diag=.true. \n / \n&test_reduction_methods_nml \n test_case = 2 \n mask_case = 2 \n \n/" | cat > input.nml
+test_expect_success "Running diag_manager with "avg" reduction method with halo output with real mask (test $my_test_count)" '
+  mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "avg" reduction method with halo output with real mask (test $my_test_count)" '
+  mpirun -n 1 ../check_time_avg
+'
+fi
+test_done