diff --git a/diag_manager/diag_data.F90 b/diag_manager/diag_data.F90
index 4e8b774afc..c601c877a9 100644
--- a/diag_manager/diag_data.F90
+++ b/diag_manager/diag_data.F90
@@ -107,7 +107,7 @@ MODULE diag_data_mod
!! to indicate to use the full axis instead of a sub-axis
INTEGER, PARAMETER :: GLO_REG_VAL_ALT = -1 !< Alternate value used in the region specification of the
!! diag_table to indicate to use the full axis instead of a sub-axis
- REAL, PARAMETER :: CMOR_MISSING_VALUE = 1.0e20 !< CMOR standard missing value
+ REAL(r8_kind), PARAMETER :: CMOR_MISSING_VALUE = 1.0e20 !< CMOR standard missing value
INTEGER, PARAMETER :: DIAG_FIELD_NOT_FOUND = -1 !< Return value for a diag_field that isn't found in the diag_table
INTEGER, PARAMETER :: latlon_gridtype = 1
INTEGER, PARAMETER :: index_gridtype = 2
@@ -390,9 +390,13 @@ MODULE diag_data_mod
!! the default behavior will do the average between day1 hour3 to day2 hour3
!
- REAL :: FILL_VALUE = NF_FILL_REAL !< Fill value used. Value will be NF90_FILL_REAL if using the
+#ifdef use_netCDF
+ REAL(r8_kind) :: FILL_VALUE = NF_FILL_REAL !< Fill value used. Value will be NF90_FILL_REAL if using the
!! netCDF module, otherwise will be 9.9692099683868690e+36.
! from file /usr/local/include/netcdf.inc
+#else
+ REAL(r8_kind) :: FILL_VALUE = 9.9692099683868690e+36
+#endif
!! @note `pack_size` and `pack_size_str` are set in diag_manager_init depending on how FMS was compiled
!! if FMS was compiled with default reals as 64bit, it will be set to 1 and "double",
@@ -405,8 +409,8 @@ MODULE diag_data_mod
!! set to "double" or "float"
!
- REAL :: EMPTY = 0.0
- REAL :: MAX_VALUE, MIN_VALUE
+ REAL(r8_kind) :: EMPTY = 0.0
+ REAL(r8_kind) :: MAX_VALUE, MIN_VALUE
!
TYPE(time_type) :: diag_init_time !< Time diag_manager_init called. If init_time not included in
diff --git a/diag_manager/diag_manager.F90 b/diag_manager/diag_manager.F90
index c1bf80dd1a..5b5357b514 100644
--- a/diag_manager/diag_manager.F90
+++ b/diag_manager/diag_manager.F90
@@ -1443,7 +1443,7 @@ LOGICAL FUNCTION send_data_0d(diag_field_id, field, time, err_msg)
TYPE(time_type), INTENT(in), OPTIONAL :: time
CHARACTER(len=*), INTENT(out), OPTIONAL :: err_msg
- REAL :: field_out(1, 1, 1) !< Local copy of field
+ CLASS(*), allocatable :: field_out(:, :, :) !< Local copy of field
! If diag_field_id is < 0 it means that this field is not registered, simply return
IF ( diag_field_id <= 0 ) THEN
@@ -1454,9 +1454,23 @@ LOGICAL FUNCTION send_data_0d(diag_field_id, field, time, err_msg)
! First copy the data to a three d array with last element 1
SELECT TYPE (field)
TYPE IS (real(kind=r4_kind))
- field_out(1, 1, 1) = field
+ allocate(real(r4_kind) :: field_out(1,1,1))
+ select type(field_out)
+ type is (real(r4_kind))
+ field_out(1, 1, 1) = field
+ class default
+ call error_mesg('diag_manager_mod::send_data_0d', &
+ & 'Error allocating field out as real(r4_kind)', FATAL)
+ end select
TYPE IS (real(kind=r8_kind))
- field_out(1, 1, 1) = real(field)
+ allocate(real(r8_kind) :: field_out(1,1,1))
+ select type(field_out)
+ type is (real(r8_kind))
+ field_out(1, 1, 1) = field
+ class default
+ call error_mesg('diag_manager_mod::send_data_0d', &
+ & 'Error allocating field out as real(r8_kind)', FATAL)
+ end select
CLASS DEFAULT
CALL error_mesg ('diag_manager_mod::send_data_0d',&
& 'The field is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
@@ -1476,7 +1490,7 @@ LOGICAL FUNCTION send_data_1d(diag_field_id, field, time, is_in, mask, rmask, ie
LOGICAL, INTENT(in), DIMENSION(:), OPTIONAL :: mask
CHARACTER(len=*), INTENT(out), OPTIONAL :: err_msg
- REAL, DIMENSION(SIZE(field(:)), 1, 1) :: field_out !< Local copy of field
+ CLASS(*), ALLOCATABLE :: field_out(:,:,:) !< Local copy of field
LOGICAL, DIMENSION(SIZE(field(:)), 1, 1) :: mask_out !< Local copy of mask
! If diag_field_id is < 0 it means that this field is not registered, simply return
@@ -1486,11 +1500,26 @@ LOGICAL FUNCTION send_data_1d(diag_field_id, field, time, is_in, mask, rmask, ie
END IF
! First copy the data to a three d array with last element 1
+ ! type checking done in diag_send_data
SELECT TYPE (field)
TYPE IS (real(kind=r4_kind))
- field_out(:, 1, 1) = field
+ allocate(real(r4_kind) :: field_out(SIZE(field),1,1))
+ select type(field_out)
+ type is (real(r4_kind))
+ field_out(:, 1, 1) = field
+ class default
+ call error_mesg('diag_manager_mod::send_data_1d', &
+ & 'Error allocating field out as real(r4_kind)', FATAL)
+ end select
TYPE IS (real(kind=r8_kind))
- field_out(:, 1, 1) = real(field)
+ allocate(real(r8_kind) :: field_out(SIZE(field),1,1))
+ select type(field_out)
+ type is (real(r8_kind))
+ field_out(:, 1, 1) = field
+ class default
+ call error_mesg('diag_manager_mod::send_data_1d', &
+ & 'Error allocating field out as real(r8_kind)', FATAL)
+ end select
CLASS DEFAULT
CALL error_mesg ('diag_manager_mod::send_data_1d',&
& 'The field is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
@@ -1545,7 +1574,7 @@ LOGICAL FUNCTION send_data_2d(diag_field_id, field, time, is_in, js_in, &
CLASS(*), INTENT(in), DIMENSION(:,:),OPTIONAL :: rmask
CHARACTER(len=*), INTENT(out), OPTIONAL :: err_msg
- REAL, DIMENSION(SIZE(field,1),SIZE(field,2),1) :: field_out !< Local copy of field
+ CLASS(*), ALLOCATABLE :: field_out(:,:,:) !< Local copy of field
LOGICAL, DIMENSION(SIZE(field,1),SIZE(field,2),1) :: mask_out !< Local copy of mask
! If diag_field_id is < 0 it means that this field is not registered, simply return
@@ -1557,9 +1586,23 @@ LOGICAL FUNCTION send_data_2d(diag_field_id, field, time, is_in, js_in, &
! First copy the data to a three d array with last element 1
SELECT TYPE (field)
TYPE IS (real(kind=r4_kind))
- field_out(:, :, 1) = field
+ allocate(real(r4_kind) :: field_out(SIZE(field,1),SIZE(field,2),1))
+ select type(field_out)
+ type is (real(r4_kind))
+ field_out(:, :, 1) = field
+ class default
+ call error_mesg('diag_manager_mod::send_data_2d', &
+ & 'Error allocating field out as real(r4_kind)', FATAL)
+ end select
TYPE IS (real(kind=r8_kind))
- field_out(:, :, 1) = real(field)
+ allocate(real(r8_kind) :: field_out(SIZE(field,1),SIZE(field,2),1))
+ select type(field_out)
+ type is (real(r8_kind))
+ field_out(:, :, 1) = field
+ class default
+ call error_mesg('diag_manager_mod::send_data_2d', &
+ & 'Error allocating field out as real(r8_kind)', FATAL)
+ end select
CLASS DEFAULT
CALL error_mesg ('diag_manager_mod::send_data_2d',&
& 'The field is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
diff --git a/diag_manager/fms_diag_object.F90 b/diag_manager/fms_diag_object.F90
index 10c8514479..894d8023d8 100644
--- a/diag_manager/fms_diag_object.F90
+++ b/diag_manager/fms_diag_object.F90
@@ -517,7 +517,7 @@ logical function fms_diag_accept_data (this, diag_field_id, field_data, mask, rm
!! the calculationslater. \note This is experimental
character(len=128) :: error_string !< Store error text
logical :: data_buffer_is_allocated !< .true. if the data buffer is allocated
- character(len=128) :: field_info !< String holding info about the field to append to the
+ character(len=256) :: field_info !< String holding info about the field to append to the
!! error message
logical, allocatable, dimension(:,:,:,:) :: oor_mask !< Out of range mask
real(kind=r8_kind) :: field_weight !< Weight to use when averaging (it will be converted
@@ -899,6 +899,11 @@ function fms_diag_do_reduction(this, field_data, diag_field_id, oor_mask, weight
return
endif
case (time_sum)
+ error_msg = buffer_ptr%do_time_sum_wrapper(field_data, oor_mask, field_ptr%get_mask_variant(), &
+ bounds_in, bounds_out, missing_value)
+ if (trim(error_msg) .ne. "") then
+ return
+ endif
case (time_average)
case (time_power)
case (time_rms)
diff --git a/diag_manager/fms_diag_output_buffer.F90 b/diag_manager/fms_diag_output_buffer.F90
index b2c54f1387..b8da2b3a62 100644
--- a/diag_manager/fms_diag_output_buffer.F90
+++ b/diag_manager/fms_diag_output_buffer.F90
@@ -34,7 +34,7 @@ module fms_diag_output_buffer_mod
use fms2_io_mod, only: FmsNetcdfFile_t, write_data, FmsNetcdfDomainFile_t, FmsNetcdfUnstructuredDomainFile_t
use fms_diag_yaml_mod, only: diag_yaml
use fms_diag_bbox_mod, only: fmsDiagIbounds_type
-use fms_diag_reduction_methods_mod, only: do_time_none, do_time_min, do_time_max
+use fms_diag_reduction_methods_mod, only: do_time_none, do_time_min, do_time_max, do_time_sum_update
use fms_diag_time_utils_mod, only: diag_time_inc
implicit none
@@ -48,10 +48,8 @@ module fms_diag_output_buffer_mod
class(*), allocatable :: buffer(:,:,:,:,:) !< 5D numeric data array
integer :: ndim !< Number of dimensions for each variable
integer, allocatable :: buffer_dims(:) !< holds the size of each dimension in the buffer
- real(r8_kind), allocatable :: counter(:,:,:,:,:) !< (x,y,z, time-of-day) used in the time averaging functions
+ real(r8_kind) :: weight_sum !< (x,y,z, time-of-day) used in the time averaging functions
integer, allocatable :: num_elements(:) !< used in time-averaging
- real(r8_kind), allocatable :: count_0d(:) !< used in time-averaging along with
- !! counter which is stored in the child types (bufferNd)
integer, allocatable :: axis_ids(:) !< Axis ids for the buffer
integer :: field_id !< The id of the field the buffer belongs to
integer :: yaml_id !< The id of the yaml id the buffer belongs to
@@ -78,6 +76,7 @@ module fms_diag_output_buffer_mod
procedure :: do_time_none_wrapper
procedure :: do_time_min_wrapper
procedure :: do_time_max_wrapper
+ procedure :: do_time_sum_wrapper
end type fmsDiagOutputBuffer_type
@@ -124,9 +123,7 @@ subroutine flush_buffer(this)
this%yaml_id = diag_null
if (allocated(this%buffer)) deallocate(this%buffer)
if (allocated(this%buffer_dims)) deallocate(this%buffer_dims)
- if (allocated(this%counter)) deallocate(this%counter)
if (allocated(this%num_elements)) deallocate(this%num_elements)
- if (allocated(this%count_0d)) deallocate(this%count_0d)
if (allocated(this%axis_ids)) deallocate(this%axis_ids)
end subroutine flush_buffer
@@ -154,38 +151,22 @@ subroutine allocate_buffer(this, buff_type, ndim, buff_sizes, field_name, diurna
type is (integer(kind=i4_kind))
allocate(integer(kind=i4_kind) :: this%buffer(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
& buff_sizes(5)))
- allocate(this%counter(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
- & buff_sizes(5)))
- allocate(this%count_0d(n_samples))
- this%counter = 0.0_r4_kind
- this%count_0d = 0.0_r4_kind
+ this%weight_sum = 0.0_r4_kind
this%buffer_type = i4
type is (integer(kind=i8_kind))
allocate(integer(kind=i8_kind) :: this%buffer(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
& buff_sizes(5)))
- allocate(this%counter(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
- & buff_sizes(5)))
- allocate(this%count_0d(n_samples))
- this%counter = 0.0_r8_kind
- this%count_0d = 0.0_r8_kind
+ this%weight_sum = 0.0_r8_kind
this%buffer_type = i8
type is (real(kind=r4_kind))
allocate(real(kind=r4_kind) :: this%buffer(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
& buff_sizes(5)))
- allocate(this%counter(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
- & buff_sizes(5)))
- allocate(this%count_0d(n_samples))
- this%counter = 0.0_r4_kind
- this%count_0d = 0.0_r4_kind
+ this%weight_sum = 0.0_r4_kind
this%buffer_type = r4
type is (real(kind=r8_kind))
allocate(real(kind=r8_kind) :: this%buffer(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
& buff_sizes(5)))
- allocate(this%counter(buff_sizes(1),buff_sizes(2),buff_sizes(3),buff_sizes(4), &
- & buff_sizes(5)))
- allocate(this%count_0d(n_samples))
- this%counter = 0.0_r8_kind
- this%count_0d = 0.0_r8_kind
+ this%weight_sum = 0.0_r8_kind
this%buffer_type = r8
class default
call mpp_error("allocate_buffer", &
@@ -194,7 +175,6 @@ subroutine allocate_buffer(this, buff_type, ndim, buff_sizes, field_name, diurna
end select
allocate(this%num_elements(n_samples))
this%num_elements = 0
- this%count_0d = 0
this%done_with_math = .false.
allocate(this%buffer_dims(5))
this%buffer_dims(1) = buff_sizes(1)
@@ -571,5 +551,42 @@ function do_time_max_wrapper(this, field_data, mask, is_masked, bounds_in, bound
end select
end select
end function do_time_max_wrapper
+
+!> @brief Does the time_sum reduction method on the buffer object
+!! @return Error message if the math was not successful
+function do_time_sum_wrapper(this, field_data, mask, is_masked, bounds_in, bounds_out, missing_value) &
+ result(err_msg)
+ class(fmsDiagOutputBuffer_type), intent(inout) :: this !< buffer object to write
+ class(*), intent(in) :: field_data(:,:,:,:) !< Buffer data for current time
+ type(fmsDiagIbounds_type), intent(in) :: bounds_in !< Indicies for the buffer passed in
+ type(fmsDiagIbounds_type), intent(in) :: bounds_out !< Indicies for the output buffer
+ logical, intent(in) :: mask(:,:,:,:) !< Mask for the field
+ logical, intent(in) :: is_masked !< .True. if the field has a mask
+ real(kind=r8_kind), intent(in) :: missing_value !< Missing_value for data points that are masked
+ character(len=50) :: err_msg
+
+ !TODO This will be expanded for integers
+ err_msg = ""
+ select type (output_buffer => this%buffer)
+ type is (real(kind=r8_kind))
+ select type (field_data)
+ type is (real(kind=r8_kind))
+ call do_time_sum_update(output_buffer, this%weight_sum, field_data, mask, is_masked, &
+ bounds_in, bounds_out, missing_value)
+ class default
+ err_msg="do_time_sum_wrapper::the output buffer and the buffer send in are not of the same type (r8_kind)"
+ end select
+ type is (real(kind=r4_kind))
+ select type (field_data)
+ type is (real(kind=r4_kind))
+ call do_time_sum_update(output_buffer, this%weight_sum, field_data, mask, is_masked, bounds_in, bounds_out, &
+ real(missing_value, kind=r4_kind))
+ class default
+ err_msg="do_time_sum_wrapper::the output buffer and the buffer send in are not of the same type (r4_kind)"
+ end select
+ class default
+ err_msg="do_time_sum_wrapper::the output buffer is not a valid type, must be real(r8_kind) or real(r4_kind)"
+ end select
+end function do_time_sum_wrapper
#endif
end module fms_diag_output_buffer_mod
diff --git a/diag_manager/fms_diag_reduction_methods.F90 b/diag_manager/fms_diag_reduction_methods.F90
index c48f9b21cd..c3d939b0f6 100644
--- a/diag_manager/fms_diag_reduction_methods.F90
+++ b/diag_manager/fms_diag_reduction_methods.F90
@@ -35,7 +35,7 @@ module fms_diag_reduction_methods_mod
private
public :: check_indices_order, init_mask, set_weight
- public :: do_time_none, do_time_min, do_time_max
+ public :: do_time_none, do_time_min, do_time_max, do_time_sum_update
!> @brief Does the time_none reduction method. See include/fms_diag_reduction_methods.inc
!TODO This needs to be extended to integers
@@ -55,6 +55,13 @@ module fms_diag_reduction_methods_mod
module procedure do_time_max_r4, do_time_max_r8
end interface do_time_max
+ !> @brief Sum update updates the buffer for any reductions that involve summation
+ !! (ie. time_sum, avg, rms, pow)
+ !!TODO This needs to be extended to integers
+ interface do_time_sum_update
+ module procedure do_time_sum_update_r4, do_time_sum_update_r8
+ end interface
+
contains
!> @brief Checks improper combinations of is, ie, js, and je.
diff --git a/diag_manager/include/fms_diag_reduction_methods.inc b/diag_manager/include/fms_diag_reduction_methods.inc
index 72332d650e..c847817724 100644
--- a/diag_manager/include/fms_diag_reduction_methods.inc
+++ b/diag_manager/include/fms_diag_reduction_methods.inc
@@ -17,6 +17,11 @@
!* License along with FMS. If not, see .
!***********************************************************************
+! for any debug prints
+#ifndef DEBUG_REDUCT
+#define DEBUG_REDUCT .true.
+#endif
+
!> @brief Do the time_none reduction method (i.e copy the correct portion of the input data)
subroutine DO_TIME_NONE_ (data_out, data_in, mask, is_masked, bounds_in, bounds_out, missing_value)
real(FMS_TRM_KIND_), intent(inout) :: data_out(:,:,:,:,:) !< output data
@@ -201,4 +206,94 @@ subroutine DO_TIME_MAX_ (data_out, data_in, mask, is_masked, bounds_in, bounds_o
enddo
enddo
endif
-end subroutine DO_TIME_MAX_
\ No newline at end of file
+end subroutine DO_TIME_MAX_
+
+!> Update the output buffer for reductions that involve summation (sum, avg, rms, pow).
+!! Elements of the running field output buffer (data_out) are set with the following:
+!!
+!! buffer(l) = buffer(l) + (weight * field(l)) ^ pow
+!!
+!! 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)
+ 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
+ logical, intent(in) :: mask(:,:,:,:) !< mask
+ logical, intent(in) :: is_masked !< .True. if the field is using a mask
+ type(fmsDiagIbounds_type), intent(in) :: bounds_in !< indices indicating the correct portion
+ !! of the input buffer
+ 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
+ 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
+
+ integer :: is_in, ie_in, js_in, je_in, ks_in, ke_in !< Starting and ending indices of each dimention for
+ !! the input buffer
+ integer :: is_out, ie_out, js_out, je_out, ks_out, ke_out !< Starting and ending indices of each dimention for
+ !! the output buffer
+ integer :: i, j, k, l !< For looping
+ real(FMS_TRM_KIND_) :: weight_loc, pow_loc !< local copies of optional arguments
+ integer, parameter :: kindl = FMS_TRM_KIND_ !< real kind size as set by macro
+
+ if(present(weight)) then
+ weight_loc = weight
+ else
+ weight_loc = 1.0_kindl
+ endif
+
+ if(present(pow)) then
+ pow_loc = weight
+ else
+ pow_loc = 1.0_kindl
+ endif
+
+ ! update with given weight for average before write
+ weight_sum = weight_sum + weight_loc
+
+ is_out = bounds_out%get_imin()
+ ie_out = bounds_out%get_imax()
+ js_out = bounds_out%get_jmin()
+ je_out = bounds_out%get_jmax()
+ ks_out = bounds_out%get_kmin()
+ ke_out = bounds_out%get_kmax()
+
+ is_in = bounds_in%get_imin()
+ ie_in = bounds_in%get_imax()
+ js_in = bounds_in%get_jmin()
+ je_in = bounds_in%get_jmax()
+ ks_in = bounds_in%get_kmin()
+ ke_in = bounds_in%get_kmax()
+
+ !> Seperated this loops for performance. If is_masked = .false. (i.e "mask" and "rmask" were never passed in)
+ !! then mask will always be .True. so the if (mask) is redudant.
+ ! TODO check if performance gain by not doing weight and pow if not needed
+ if (is_masked) then
+ do k = 0, ke_out - ks_out
+ do j = 0, je_out - js_out
+ do i = 0, ie_out - is_out
+ where (mask(is_in + i, js_in + j, ks_in + k, :))
+ 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
+ elsewhere
+ data_out(is_out + i, js_out + j, ks_out + k, :, 1) = missing_value
+ endwhere
+ enddo
+ enddo
+ enddo
+ else
+ ! doesn't need to loop through l if no mask, just sums the 1d slices
+ 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_in(is_in +i, js_in + j, ks_in + k, :) * weight_loc) ** pow_loc
+ enddo
+ enddo
+ enddo
+ endif
+end subroutine DO_TIME_SUM_UPDATE_
diff --git a/diag_manager/include/fms_diag_reduction_methods_r4.fh b/diag_manager/include/fms_diag_reduction_methods_r4.fh
index c3bc29296a..a3c499b12e 100644
--- a/diag_manager/include/fms_diag_reduction_methods_r4.fh
+++ b/diag_manager/include/fms_diag_reduction_methods_r4.fh
@@ -35,6 +35,9 @@
#undef DO_TIME_MAX_
#define DO_TIME_MAX_ do_time_max_r4
+#undef DO_TIME_SUM_UPDATE_
+#define DO_TIME_SUM_UPDATE_ do_time_sum_update_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 a3e3d68376..d550293113 100644
--- a/diag_manager/include/fms_diag_reduction_methods_r8.fh
+++ b/diag_manager/include/fms_diag_reduction_methods_r8.fh
@@ -35,6 +35,9 @@
#undef DO_TIME_MAX_
#define DO_TIME_MAX_ do_time_max_r8
+#undef DO_TIME_SUM_UPDATE_
+#define DO_TIME_SUM_UPDATE_ do_time_sum_update_r8
+
#include "fms_diag_reduction_methods.inc"
!> @}
diff --git a/test_fms/diag_manager/Makefile.am b/test_fms/diag_manager/Makefile.am
index de682cc7ee..35c0aa3198 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_min check_time_max check_time_sum
# This is the source code for the test.
test_diag_manager_SOURCES = test_diag_manager.F90
@@ -47,18 +47,20 @@ test_reduction_methods_SOURCES = testing_utils.F90 test_reduction_methods.F90
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
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
+TESTS = test_diag_manager2.sh test_time_none.sh test_time_min.sh test_time_max.sh test_time_sum.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
+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
if USING_YAML
skipflag=""
diff --git a/test_fms/diag_manager/check_time_sum.F90 b/test_fms/diag_manager/check_time_sum.F90
new file mode 100644
index 0000000000..03d38f21a2
--- /dev/null
+++ b/test_fms/diag_manager/check_time_sum.F90
@@ -0,0 +1,264 @@
+!***********************************************************************
+!* 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_sum" reduction method
+program check_time_sum
+ 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_sum.nc", "read")) &
+ call mpp_error(FATAL, "unable to open test_sum.nc")
+
+ if (.not. open_file(fileobj1, "test_sum_regional.nc.0004", "read")) &
+ call mpp_error(FATAL, "unable to open test_sum_regional.nc.0004")
+
+ if (.not. open_file(fileobj2, "test_sum_regional.nc.0005", "read")) &
+ call mpp_error(FATAL, "unable to open test_sum_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_sum - time_level:", string(ti)
+ call read_data(fileobj, "var0_sum", 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_sum - time_level:", string(ti)
+ call read_data(fileobj, "var1_sum", 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_sum - time_level:", string(ti)
+ call read_data(fileobj, "var2_sum", 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_sum - time_level:", string(ti)
+ call read_data(fileobj, "var3_sum", 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 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_sum in the first regional file- time_level:", string(ti)
+ call read_data(fileobj1, "var3_sum", 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_sum in the second regional file- time_level:", string(ti)
+ call read_data(fileobj2, "var3_sum", 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_sum = 0 !< sum of time step increments to use in generating reference data
+
+ ! sums integers for decimal part of field input
+ ! ie. level 1 = 1+2+..+6
+ ! 2 = 7+8+..+12
+ step_sum = 0
+ do i=(time_level-1)*file_freq+1, time_level*file_freq
+ step_sum = step_sum + i
+ enddo
+
+ ! 0d answer is:
+ ! (1011 * frequency sum'd over )
+ ! + ( 1/100 * sum of time step increments )
+ buffer_exp = real((1000.0_r8_kind+10.0_r8_kind+1.0_r8_kind) * file_freq + &
+ real(step_sum,r8_kind)/100.0_r8_kind, kind=r4_kind)
+
+ if (abs(buffer - buffer_exp) > 0.0) then
+ print *, mpp_pe(), time_level, buffer_exp, buffer
+ call mpp_error(FATAL, "Check_time_sum::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 !< sum 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)
+ do ii = 1, size(buffer, 1)
+ buffer_exp = 0.0
+ ! fails with both precisions
+ !do n=(time_level-1)*file_freq+1, time_level*file_freq
+ ! buffer_exp = real(buffer_exp + 1000.0_r8_kind * ii + 11.0_r8_kind + (n/100.0_r8_kind), r4_kind)
+ !enddo
+ ! passes with r8 defaults, fails with r4
+ buffer_exp = real( &
+ file_freq * (real(ii, kind=r8_kind)*1000.0_r8_kind +10.0_r8_kind+1.0_r8_kind) + &
+ real(step_sum, kind=r8_kind)/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, "sum of time steps:", step_sum
+ print *, "diff:", abs(buffer(ii) - buffer_exp)
+ call mpp_error(FATAL, "Check_time_sum::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_sum !< sum of time step increments to use in generating reference data
+
+ step_sum = 0
+ do i=(time_level-1)*file_freq+1, time_level*file_freq
+ step_sum = step_sum + i
+ enddo
+
+ ! 2d answer is
+ ! ((i * 1000 + j * 10 + 1) * frequency) + (sum of time steps)
+ do ii = 1, size(buffer, 1)
+ do j = 1, size(buffer, 2)
+ buffer_exp = real(real(ii, kind=r8_kind)* 6000.0_kindl+ &
+ 60.0_kindl*real(j, kind=r8_kind)+6.0_kindl + &
+ real(step_sum, kind=r8_kind)/100_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 *, mpp_pe(), ii, j, buffer(ii, j), buffer_exp
+ call mpp_error(FATAL, "Check_time_sum::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_sum!< sum of time step increments to use in generating reference data
+
+ step_sum = 0
+ do i=(time_level-1)*file_freq+1, time_level*file_freq
+ step_sum = step_sum + 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) + (sum 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)* 6000.0_kindl + &
+ 60.0_kindl*real(j+ny_oset, kind=r8_kind) + &
+ 6.0_kindl*real(k+nz_oset, kind=r8_kind) + &
+ real(step_sum, kind=r8_kind)/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(), ii, j, k, buffer(ii, j, k), buffer_exp
+ call mpp_error(FATAL, "Check_time_sum::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_time_sum.sh b/test_fms/diag_manager/test_time_sum.sh
new file mode 100755
index 0000000000..18f923cbb4
--- /dev/null
+++ b/test_fms/diag_manager/test_time_sum.sh
@@ -0,0 +1,166 @@
+#!/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_sum
+base_date: 2 1 1 0 0 0
+diag_files:
+- file_name: test_sum
+ time_units: hours
+ unlimdim: time
+ freq: 6 hours
+ varlist:
+ - module: ocn_mod
+ var_name: var0
+ output_name: var0_sum
+ reduction: sum
+ kind: r4
+ - module: ocn_mod
+ var_name: var1
+ output_name: var1_sum
+ reduction: sum
+ kind: r4
+ - module: ocn_mod
+ var_name: var2
+ output_name: var2_sum
+ reduction: sum
+ kind: r4
+ - module: ocn_mod
+ var_name: var3
+ output_name: var3_sum
+ reduction: sum
+ kind: r4
+ - module: ocn_mod
+ var_name: var3
+ output_name: var3_Z
+ reduction: sum
+ zbounds: 2. 3.
+ kind: r4
+- file_name: test_sum_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_sum
+ reduction: sum
+ zbounds: 2. 3.
+ kind: r4
+_EOF
+
+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 "sum" reduction method (test $my_test_count)" '
+ mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "sum" reduction method (test $my_test_count)" '
+ mpirun -n 1 ../check_time_sum
+'
+
+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 "sum" reduction method, logical mask (test $my_test_count)" '
+ mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "sum" reduction method, logical mask (test $my_test_count)" '
+ mpirun -n 1 ../check_time_sum
+'
+
+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 "sum" reduction method, real mask (test $my_test_count)" '
+ mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "sum" reduction method, real mask (test $my_test_count)" '
+ mpirun -n 1 ../check_time_sum
+'
+
+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 = 1 \n \n/" | cat > input.nml
+test_expect_success "Running diag_manager with "sum" reduction method with openmp (test $my_test_count)" '
+ mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "sum" reduction method with openmp (test $my_test_count)" '
+ mpirun -n 1 ../check_time_sum
+'
+
+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 "sum" reduction method with openmp, logical mask (test $my_test_count)" '
+ mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "sum" reduction method with openmp, logical mask (test $my_test_count)" '
+ mpirun -n 1 ../check_time_sum
+'
+
+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 "sum" reduction method with openmp, real mask (test $my_test_count)" '
+ mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "sum" reduction method with openmp, real mask (test $my_test_count)" '
+ mpirun -n 1 ../check_time_sum
+'
+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 = 2 \n \n/" | cat > input.nml
+test_expect_success "Running diag_manager with "sum" reduction method with halo output (test $my_test_count)" '
+ mpirun -n 6 ../test_reduction_methods
+'
+test_expect_success "Checking answers for the "sum" reduction method with halo output (test $my_test_count)" '
+ mpirun -n 1 ../check_time_sum
+'
+
+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 "sum" 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 "sum" reduction method with halo output with logical mask (test $my_test_count)" '
+ mpirun -n 1 ../check_time_sum
+'
+
+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 "sum" 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 "sum" reduction method with halo output with real mask (test $my_test_count)" '
+ mpirun -n 1 ../check_time_sum
+'
+fi
+test_done