dmUpdate: add time_average reduction method and test (NOAA-GFDL#1421)

* add tests and routines from other branch

* move weight check to output buffer, more comments + cleanup

* simplify reduction done

* update avg test for send_data_r4

* Add time to the buffer object so the counter is only increased if data has not been recived for that time

---------

Co-authored-by: rem1776 <[email protected]>
Co-authored-by: Uriel Ramirez <[email protected]>

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)