From 9048fa47315219f99fe5c61590444a0996dd5eff Mon Sep 17 00:00:00 2001
From: Ryan Mulhall <35538242+rem1776@users.noreply.github.com>
Date: Wed, 27 Dec 2023 11:25:38 -0500
Subject: [PATCH] feat: modern diag add time_average reduction method and test
(#1421)
---
diag_manager/diag_data.F90 | 12 +-
diag_manager/fms_diag_axis_object.F90 | 5 +-
diag_manager/fms_diag_field_object.F90 | 53 +++-
diag_manager/fms_diag_file_object.F90 | 18 ++
diag_manager/fms_diag_object.F90 | 74 ++++-
diag_manager/fms_diag_output_buffer.F90 | 84 +++++-
diag_manager/fms_diag_reduction_methods.F90 | 10 +-
.../include/fms_diag_reduction_methods.inc | 37 ++-
.../include/fms_diag_reduction_methods_r4.fh | 3 +
.../include/fms_diag_reduction_methods_r8.fh | 3 +
test_fms/diag_manager/Makefile.am | 8 +-
test_fms/diag_manager/check_time_avg.F90 | 270 ++++++++++++++++++
test_fms/diag_manager/test_diag_manager2.sh | 2 +-
test_fms/diag_manager/test_modern_diag.F90 | 2 +-
test_fms/diag_manager/test_time_avg.sh | 180 ++++++++++++
15 files changed, 725 insertions(+), 36 deletions(-)
create mode 100644 test_fms/diag_manager/check_time_avg.F90
create mode 100755 test_fms/diag_manager/test_time_avg.sh
diff --git a/diag_manager/diag_data.F90 b/diag_manager/diag_data.F90
index 56aa041fdc..f73776ff75 100644
--- a/diag_manager/diag_data.F90
+++ b/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
diff --git a/diag_manager/fms_diag_axis_object.F90 b/diag_manager/fms_diag_axis_object.F90
index 8f22f7d2db..e74ccabeff 100644
--- a/diag_manager/fms_diag_axis_object.F90
+++ b/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)
diff --git a/diag_manager/fms_diag_field_object.F90 b/diag_manager/fms_diag_field_object.F90
index e723ce8410..65fd44719b 100644
--- a/diag_manager/fms_diag_field_object.F90
+++ b/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
diff --git a/diag_manager/fms_diag_file_object.F90 b/diag_manager/fms_diag_file_object.F90
index 687f609252..1aa5baf899 100644
--- a/diag_manager/fms_diag_file_object.F90
+++ b/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
diff --git a/diag_manager/fms_diag_object.F90 b/diag_manager/fms_diag_object.F90
index d92d6a9cf2..0d7cddf37e 100644
--- a/diag_manager/fms_diag_object.F90
+++ b/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,6 +708,8 @@ 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,
@@ -706,11 +717,21 @@ subroutine fms_diag_do_io(this, is_end_of_run)
#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)
diff --git a/diag_manager/fms_diag_output_buffer.F90 b/diag_manager/fms_diag_output_buffer.F90
index b8da2b3a62..eed366bee1 100644
--- a/diag_manager/fms_diag_output_buffer.F90
+++ b/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,7 +608,7 @@ 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
@@ -580,7 +616,7 @@ function do_time_sum_wrapper(this, field_data, mask, is_masked, bounds_in, bound
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
@@ -588,5 +624,39 @@ function do_time_sum_wrapper(this, field_data, mask, is_masked, bounds_in, bound
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
diff --git a/diag_manager/fms_diag_reduction_methods.F90 b/diag_manager/fms_diag_reduction_methods.F90
index c3d939b0f6..801f6ba557 100644
--- a/diag_manager/fms_diag_reduction_methods.F90
+++ b/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.
diff --git a/diag_manager/include/fms_diag_reduction_methods.inc b/diag_manager/include/fms_diag_reduction_methods.inc
index c847817724..c443a945b4 100644
--- a/diag_manager/include/fms_diag_reduction_methods.inc
+++ b/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,7 +291,7 @@ 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
@@ -297,3 +299,30 @@ subroutine DO_TIME_SUM_UPDATE_(data_out, weight_sum, data_in, mask, is_masked, b
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
+
diff --git a/diag_manager/include/fms_diag_reduction_methods_r4.fh b/diag_manager/include/fms_diag_reduction_methods_r4.fh
index a3c499b12e..04a4f4f0ba 100644
--- a/diag_manager/include/fms_diag_reduction_methods_r4.fh
+++ b/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"
!> @}
diff --git a/diag_manager/include/fms_diag_reduction_methods_r8.fh b/diag_manager/include/fms_diag_reduction_methods_r8.fh
index d550293113..bff7f44ac2 100644
--- a/diag_manager/include/fms_diag_reduction_methods_r8.fh
+++ b/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"
!> @}
diff --git a/test_fms/diag_manager/Makefile.am b/test_fms/diag_manager/Makefile.am
index 35c0aa3198..ea251f291f 100644
--- a/test_fms/diag_manager/Makefile.am
+++ b/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=""
diff --git a/test_fms/diag_manager/check_time_avg.F90 b/test_fms/diag_manager/check_time_avg.F90
new file mode 100644
index 0000000000..6a1d527537
--- /dev/null
+++ b/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 .
+!***********************************************************************
+
+!> @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
diff --git a/test_fms/diag_manager/test_diag_manager2.sh b/test_fms/diag_manager/test_diag_manager2.sh
index 813e225156..449c3ed82f 100755
--- a/test_fms/diag_manager/test_diag_manager2.sh
+++ b/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
diff --git a/test_fms/diag_manager/test_modern_diag.F90 b/test_fms/diag_manager/test_modern_diag.F90
index 8205b8eee1..b39eb44594 100644
--- a/test_fms/diag_manager/test_modern_diag.F90
+++ b/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
diff --git a/test_fms/diag_manager/test_time_avg.sh b/test_fms/diag_manager/test_time_avg.sh
new file mode 100755
index 0000000000..7c9752231c
--- /dev/null
+++ b/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 .
+#***********************************************************************
+
+# 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