! CSIRO Open Source Software License Agreement (variation of the BSD / MIT License) ! Copyright (c) 2015, Commonwealth Scientific and Industrial Research Organisation ! (CSIRO) ABN 41 687 119 230. ! TODO(Sean): The preprocessor define ENFORCE_SINGLE_PRECISION is enabled ! temporarily to restore bitwise reproducibility with the previous output module ! which enforces double precision variables to be sampled using single precision ! arrays. #define ENFORCE_SINGLE_PRECISION module aggregator_mod !* This module defines the `aggregator_t` type and its extensions, which are ! used to perform various types of time aggregations (e.g., mean, sum, min, max) ! on source data arrays. ! ! Aggregators are utilised by first initialising an aggregator via ! [[new_aggregator]] which accepts as its argument a source data array to be ! sampled. The source data array argument is then associated with a pointer ! which is used to sample the source data array to compute the time aggregated ! value. The source data array argument must be an allocated pointer or ! allocatable array. The source data array argument is typically a CABLE working ! variable whos memory remains allocated throughout the duration of the ! simulation. Aggregator instances can be represented using a polymorphic ! `aggregator_t` type or a concrete type corresponding to the data type and rank ! of the source data array (e.g., `aggregator_real32_2d_t` for a 2D 32-bit real ! source data array). For example: ! ! ```fortran ! real(kind=real32), dimension(:,:), allocatable :: source_data_array ! class(aggregator_t), allocatable :: aggregator_polymorphic ! type(aggregator_real32_2d_t) :: aggregator_concrete ! ! allocate(source_data_array(42, 42)) ! ! aggregator_polymorphic = new_aggregator(source_data_array) ! aggregator_concrete = new_aggregator(source_data_array) ! ``` ! ! Aggregator instances must then be initialised by calling their `init` ! method, which will allocate the necessary memory for the aggregated data and ! set the accumulation and reset methods according to the specified aggregation ! method (e.g., mean, sum, min, max). For example: ! ! ```fortran ! call aggregator_polymorphic%init('mean') ! call aggregator_concrete%init('sum') ! ``` ! ! Once initialised, the aggregator can be accumulated any number of times ! throughout a simulation by calling its `accumulate` method, which will update ! the aggregated data from the source data using the appropriate aggregation ! method. Once the aggregator has accumulated enough values over the time ! interval of interest, the aggregated data can then be accessed via the ! `aggregated_data` component of the specific aggregator type (e.g., ! `aggregator_real32_2d_t%aggregated_data`). To reset the aggregator back to its ! initial state for the next time interval, the `reset` method can be called, ! which will reset the aggregator according to the specified aggregation method. ! ! Resources allocated by the aggregator on initialisation will automatically be ! deallocated when the aggregator instance goes out of scope. use iso_fortran_env, only: int32, real32, real64 use cable_error_handler_mod, only: cable_abort implicit none private public :: aggregator_t public :: aggregator_int32_0d_t public :: aggregator_int32_1d_t public :: aggregator_int32_2d_t public :: aggregator_int32_3d_t public :: aggregator_real32_0d_t public :: aggregator_real32_1d_t public :: aggregator_real32_2d_t public :: aggregator_real32_3d_t public :: aggregator_real64_0d_t public :: aggregator_real64_1d_t public :: aggregator_real64_2d_t public :: aggregator_real64_3d_t public :: new_aggregator type, abstract :: aggregator_t !* The `aggregator_t` type is an abstract base type for performing time ! aggregations on source data arrays. Support for different data types and ! array ranks is provided through extensions of this type. integer :: counter = 0 !! The number of times the source data has been accumulated. procedure(accumulate_data), pointer :: accumulate !* A procedure pointer to the accumulation method used to compute the time ! aggregated quantity. The specific method (e.g., mean, sum, min, max) is ! set by the `set_method` procedure. procedure(reset_data), pointer :: reset !* A procedure pointer to the reset method used to reset the aggregated data. The ! specific method (e.g., mean, sum, min, max) is set by the `set_method` procedure. contains procedure :: init => aggregator_init !! Initialise the aggregator. procedure :: type => aggregator_type !! Return a string identifier of the aggregator type. procedure :: rank => aggregator_rank !! Return the rank of the aggregator. procedure :: shape => aggregator_shape !! Return the shape of the aggregator. procedure :: scale => aggregator_scale !! Scale the aggregated data by a specified factor. procedure :: div => aggregator_div !! Divide the aggregated data by a specified factor. procedure :: offset => aggregator_offset !! Add a specified offset to the aggregated data. procedure, private :: set_method => aggregator_set_method !! Set the aggregation method. end type aggregator_t abstract interface !* Interfaces for the procedure pointers in the `aggregator_t` type to be ! implemented by the specific aggregation methods (e.g., mean, sum, min, max). subroutine accumulate_data(this, scale, div, offset) !! Accumulate the aggregated data from the source data. import aggregator_t class(aggregator_t), intent(inout) :: this real, intent(in), optional :: scale !* An optional scaling factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: div !* An optional division factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: offset !* An optional offset to add to the source data before accumulation. ! Defaults to 0.0 if not provided. end subroutine accumulate_data subroutine reset_data(this) !! Reset the aggregated data to its initial state. import aggregator_t class(aggregator_t), intent(inout) :: this end subroutine reset_data end interface type, extends(aggregator_t) :: aggregator_int32_0d_t !! An aggregator for 0-dimensional (scalar) 32-bit integer data. integer(kind=int32), allocatable :: aggregated_data integer(kind=int32), pointer :: source_data => null() end type aggregator_int32_0d_t type, extends(aggregator_t) :: aggregator_int32_1d_t !! An aggregator for 1-dimensional 32-bit integer data. integer(kind=int32), dimension(:), allocatable :: aggregated_data integer(kind=int32), dimension(:), pointer :: source_data => null() end type aggregator_int32_1d_t type, extends(aggregator_t) :: aggregator_int32_2d_t !! An aggregator for 2-dimensional 32-bit integer data. integer(kind=int32), dimension(:,:), allocatable :: aggregated_data integer(kind=int32), dimension(:,:), pointer :: source_data => null() end type aggregator_int32_2d_t type, extends(aggregator_t) :: aggregator_int32_3d_t !! An aggregator for 3-dimensional 32-bit integer data. integer(kind=int32), dimension(:,:,:), allocatable :: aggregated_data integer(kind=int32), dimension(:,:,:), pointer :: source_data => null() end type aggregator_int32_3d_t type, extends(aggregator_t) :: aggregator_real32_0d_t !! An aggregator for 0-dimensional (scalar) 32-bit real data. real(kind=real32), allocatable :: aggregated_data real(kind=real32), pointer :: source_data => null() end type aggregator_real32_0d_t type, extends(aggregator_t) :: aggregator_real32_1d_t !! An aggregator for 1-dimensional 32-bit real data. real(kind=real32), dimension(:), allocatable :: aggregated_data real(kind=real32), dimension(:), pointer :: source_data => null() end type aggregator_real32_1d_t type, extends(aggregator_t) :: aggregator_real32_2d_t !! An aggregator for 2-dimensional 32-bit real data. real(kind=real32), dimension(:,:), allocatable :: aggregated_data real(kind=real32), dimension(:,:), pointer :: source_data => null() end type aggregator_real32_2d_t type, extends(aggregator_t) :: aggregator_real32_3d_t !! An aggregator for 3-dimensional 32-bit real data. real(kind=real32), dimension(:,:,:), allocatable :: aggregated_data real(kind=real32), dimension(:,:,:), pointer :: source_data => null() end type aggregator_real32_3d_t type, extends(aggregator_t) :: aggregator_real64_0d_t !! An aggregator for 0-dimensional (scalar) 64-bit real data. #ifdef ENFORCE_SINGLE_PRECISION real(kind=real32), allocatable :: aggregated_data #else real(kind=real64), allocatable :: aggregated_data #endif real(kind=real64), pointer :: source_data => null() end type aggregator_real64_0d_t type, extends(aggregator_t) :: aggregator_real64_1d_t !! An aggregator for 1-dimensional 64-bit real data. #ifdef ENFORCE_SINGLE_PRECISION real(kind=real32), dimension(:), allocatable :: aggregated_data #else real(kind=real64), dimension(:), allocatable :: aggregated_data #endif real(kind=real64), dimension(:), pointer :: source_data => null() end type aggregator_real64_1d_t type, extends(aggregator_t) :: aggregator_real64_2d_t !! An aggregator for 2-dimensional 64-bit real data. #ifdef ENFORCE_SINGLE_PRECISION real(kind=real32), dimension(:,:), allocatable :: aggregated_data #else real(kind=real64), dimension(:,:), allocatable :: aggregated_data #endif real(kind=real64), dimension(:,:), pointer :: source_data => null() end type aggregator_real64_2d_t type, extends(aggregator_t) :: aggregator_real64_3d_t !! An aggregator for 3-dimensional 64-bit real data. #ifdef ENFORCE_SINGLE_PRECISION real(kind=real32), dimension(:,:,:), allocatable :: aggregated_data #else real(kind=real64), dimension(:,:,:), allocatable :: aggregated_data #endif real(kind=real64), dimension(:,:,:), pointer :: source_data => null() end type aggregator_real64_3d_t interface new_aggregator !* Factory interface for creating new aggregator instances. The specific ! type of aggregator created is determined by the type of the source data ! array provided. module procedure new_aggregator_int32_0d_t module procedure new_aggregator_int32_1d_t module procedure new_aggregator_int32_2d_t module procedure new_aggregator_int32_3d_t module procedure new_aggregator_real32_0d module procedure new_aggregator_real32_1d module procedure new_aggregator_real32_2d module procedure new_aggregator_real32_3d module procedure new_aggregator_real64_0d module procedure new_aggregator_real64_1d module procedure new_aggregator_real64_2d module procedure new_aggregator_real64_3d end interface contains subroutine aggregator_init(this, method) !* Initialise the aggregator by allocating the aggregated data array and its ! aggregation method. The values in the aggregated data array are reset ! according to the specified aggregation method. class(aggregator_t), intent(inout) :: this character(len=*), intent(in) :: method !! The aggregation method to use (e.g., "mean", "sum", "point", "min", "max"). select type (this) type is (aggregator_int32_0d_t) if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) type is (aggregator_int32_1d_t) if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) type is (aggregator_int32_2d_t) if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) type is (aggregator_int32_3d_t) if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) type is (aggregator_real32_0d_t) if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) type is (aggregator_real32_1d_t) if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) type is (aggregator_real32_2d_t) if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) type is (aggregator_real32_3d_t) if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) type is (aggregator_real64_0d_t) if (.not. allocated(this%aggregated_data)) allocate( & #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data, mold=real(this%source_data, kind=real32) & #else this%aggregated_data, mold=real(this%source_data, kind=real64) & #endif ) type is (aggregator_real64_1d_t) if (.not. allocated(this%aggregated_data)) allocate( & #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data, mold=real(this%source_data, kind=real32) & #else this%aggregated_data, mold=real(this%source_data, kind=real64) & #endif ) type is (aggregator_real64_2d_t) if (.not. allocated(this%aggregated_data)) allocate( & #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data, mold=real(this%source_data, kind=real32) & #else this%aggregated_data, mold=real(this%source_data, kind=real64) & #endif ) type is (aggregator_real64_3d_t) if (.not. allocated(this%aggregated_data)) allocate( & #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data, mold=real(this%source_data, kind=real32) & #else this%aggregated_data, mold=real(this%source_data, kind=real64) & #endif ) class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select call this%set_method(method) call this%reset() end subroutine aggregator_init subroutine aggregator_set_method(this, method) !* Set the aggregation method for the aggregator by assigning the appropriate ! accumulation and reset procedures based on the specified method. class(aggregator_t), intent(inout) :: this character(len=*), intent(in) :: method !! The aggregation method to use (e.g., "mean", "sum", "point", "min", "max"). if (method == "mean") then this%accumulate => mean_accumulate this%reset => other_reset elseif (method == "sum") then this%accumulate => sum_accumulate this%reset => other_reset elseif (method == "point") then this%accumulate => point_accumulate this%reset => point_reset elseif (method == "min") then this%accumulate => min_accumulate this%reset => min_reset elseif (method == "max") then this%accumulate => max_accumulate this%reset => max_reset else call cable_abort("Aggregation method "//method//" is invalid.", file=__FILE__, line=__LINE__) endif end subroutine aggregator_set_method character(16) function aggregator_type(this) !! Return a string identifier of the aggregator type (e.g., "int32", "real32", "real64"). class(aggregator_t), intent(in) :: this select type (this) type is (aggregator_int32_0d_t) aggregator_type = "int32" type is (aggregator_int32_1d_t) aggregator_type = "int32" type is (aggregator_int32_2d_t) aggregator_type = "int32" type is (aggregator_int32_3d_t) aggregator_type = "int32" type is (aggregator_real32_0d_t) aggregator_type = "real32" type is (aggregator_real32_1d_t) aggregator_type = "real32" type is (aggregator_real32_2d_t) aggregator_type = "real32" type is (aggregator_real32_3d_t) aggregator_type = "real32" type is (aggregator_real64_0d_t) aggregator_type = "real64" type is (aggregator_real64_1d_t) aggregator_type = "real64" type is (aggregator_real64_2d_t) aggregator_type = "real64" type is (aggregator_real64_3d_t) aggregator_type = "real64" class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select end function aggregator_type integer function aggregator_rank(this) !! Return the rank of the aggregator. class(aggregator_t), intent(in) :: this select type (this) type is (aggregator_int32_0d_t) aggregator_rank = 0 type is (aggregator_int32_1d_t) aggregator_rank = 1 type is (aggregator_int32_2d_t) aggregator_rank = 2 type is (aggregator_int32_3d_t) aggregator_rank = 3 type is (aggregator_real32_0d_t) aggregator_rank = 0 type is (aggregator_real32_1d_t) aggregator_rank = 1 type is (aggregator_real32_2d_t) aggregator_rank = 2 type is (aggregator_real32_3d_t) aggregator_rank = 3 type is (aggregator_real64_0d_t) aggregator_rank = 0 type is (aggregator_real64_1d_t) aggregator_rank = 1 type is (aggregator_real64_2d_t) aggregator_rank = 2 type is (aggregator_real64_3d_t) aggregator_rank = 3 class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select end function aggregator_rank function aggregator_shape(this) result(agg_shape) !! Return the shape of the aggregator. class(aggregator_t), intent(in) :: this integer, allocatable :: agg_shape(:) select type (this) type is (aggregator_int32_0d_t) agg_shape = shape(this%source_data) type is (aggregator_int32_1d_t) agg_shape = shape(this%source_data) type is (aggregator_int32_2d_t) agg_shape = shape(this%source_data) type is (aggregator_int32_3d_t) agg_shape = shape(this%source_data) type is (aggregator_real32_0d_t) agg_shape = shape(this%source_data) type is (aggregator_real32_1d_t) agg_shape = shape(this%source_data) type is (aggregator_real32_2d_t) agg_shape = shape(this%source_data) type is (aggregator_real32_3d_t) agg_shape = shape(this%source_data) type is (aggregator_real64_0d_t) agg_shape = shape(this%source_data) type is (aggregator_real64_1d_t) agg_shape = shape(this%source_data) type is (aggregator_real64_2d_t) agg_shape = shape(this%source_data) type is (aggregator_real64_3d_t) agg_shape = shape(this%source_data) class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select end function aggregator_shape subroutine aggregator_scale(this, scale) !! Scale the aggregated data by a specified factor. class(aggregator_t), intent(inout) :: this real, intent(in) :: scale !! The factor by which to scale the aggregated data. select type (this) type is (aggregator_int32_0d_t) this%aggregated_data = this%aggregated_data * scale type is (aggregator_int32_1d_t) this%aggregated_data = this%aggregated_data * scale type is (aggregator_int32_2d_t) this%aggregated_data = this%aggregated_data * scale type is (aggregator_int32_3d_t) this%aggregated_data = this%aggregated_data * scale type is (aggregator_real32_0d_t) this%aggregated_data = this%aggregated_data * scale type is (aggregator_real32_1d_t) this%aggregated_data = this%aggregated_data * scale type is (aggregator_real32_2d_t) this%aggregated_data = this%aggregated_data * scale type is (aggregator_real32_3d_t) this%aggregated_data = this%aggregated_data * scale type is (aggregator_real64_0d_t) this%aggregated_data = this%aggregated_data * scale type is (aggregator_real64_1d_t) this%aggregated_data = this%aggregated_data * scale type is (aggregator_real64_2d_t) this%aggregated_data = this%aggregated_data * scale type is (aggregator_real64_3d_t) this%aggregated_data = this%aggregated_data * scale class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select end subroutine aggregator_scale subroutine aggregator_div(this, div) !! Divide the aggregated data by a specified factor. class(aggregator_t), intent(inout) :: this real, intent(in) :: div !! The factor by which to divide the aggregated data. select type (this) type is (aggregator_int32_0d_t) this%aggregated_data = this%aggregated_data / div type is (aggregator_int32_1d_t) this%aggregated_data = this%aggregated_data / div type is (aggregator_int32_2d_t) this%aggregated_data = this%aggregated_data / div type is (aggregator_int32_3d_t) this%aggregated_data = this%aggregated_data / div type is (aggregator_real32_0d_t) this%aggregated_data = this%aggregated_data / div type is (aggregator_real32_1d_t) this%aggregated_data = this%aggregated_data / div type is (aggregator_real32_2d_t) this%aggregated_data = this%aggregated_data / div type is (aggregator_real32_3d_t) this%aggregated_data = this%aggregated_data / div type is (aggregator_real64_0d_t) this%aggregated_data = this%aggregated_data / div type is (aggregator_real64_1d_t) this%aggregated_data = this%aggregated_data / div type is (aggregator_real64_2d_t) this%aggregated_data = this%aggregated_data / div type is (aggregator_real64_3d_t) this%aggregated_data = this%aggregated_data / div class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select end subroutine aggregator_div subroutine aggregator_offset(this, offset) !! Offset the aggregated data by a specified value. class(aggregator_t), intent(inout) :: this real, intent(in) :: offset !! The value by which to offset the aggregated data. select type (this) type is (aggregator_int32_0d_t) this%aggregated_data = this%aggregated_data + offset type is (aggregator_int32_1d_t) this%aggregated_data = this%aggregated_data + offset type is (aggregator_int32_2d_t) this%aggregated_data = this%aggregated_data + offset type is (aggregator_int32_3d_t) this%aggregated_data = this%aggregated_data + offset type is (aggregator_real32_0d_t) this%aggregated_data = this%aggregated_data + offset type is (aggregator_real32_1d_t) this%aggregated_data = this%aggregated_data + offset type is (aggregator_real32_2d_t) this%aggregated_data = this%aggregated_data + offset type is (aggregator_real32_3d_t) this%aggregated_data = this%aggregated_data + offset type is (aggregator_real64_0d_t) this%aggregated_data = this%aggregated_data + offset type is (aggregator_real64_1d_t) this%aggregated_data = this%aggregated_data + offset type is (aggregator_real64_2d_t) this%aggregated_data = this%aggregated_data + offset type is (aggregator_real64_3d_t) this%aggregated_data = this%aggregated_data + offset class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select end subroutine aggregator_offset subroutine get_accumulate_args(scale, div, offset, scale_out, div_out, offset_out) !* Helper subroutine to get initialise optional scale, div, and offset arguments for ! accumulate procedures. real, intent(in), optional :: scale !* An optional scaling factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: div !* An optional division factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: offset !* An optional offset to add to the source data before accumulation. ! Defaults to 0.0 if not provided. real :: scale_out, div_out, offset_out if (present(scale)) then scale_out = scale else scale_out = 1.0 end if if (present(div)) then div_out = div else div_out = 1.0 end if if (present(offset)) then offset_out = offset else offset_out = 0.0 end if end subroutine subroutine mean_accumulate(this, scale, div, offset) !* Accumulate the aggregated data from the source data using the mean ! aggregation method. class(aggregator_t), intent(inout) :: this real, intent(in), optional :: scale !* An optional scaling factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: div !* An optional division factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: offset !* An optional offset to add to the source data before accumulation. ! Defaults to 0.0 if not provided. real :: scale_val, div_val, offset_val call get_accumulate_args(scale, div, offset, scale_val, div_val, offset_val) select type (this) type is (aggregator_real32_0d_t) this%aggregated_data = this%aggregated_data + ( & scale_val * this%source_data / div + offset_val - this%aggregated_data & ) / (this%counter + 1) type is (aggregator_real32_1d_t) this%aggregated_data = this%aggregated_data + ( & scale_val * this%source_data / div + offset_val - this%aggregated_data & ) / (this%counter + 1) type is (aggregator_real32_2d_t) this%aggregated_data = this%aggregated_data + ( & scale_val * this%source_data / div + offset_val - this%aggregated_data & ) / (this%counter + 1) type is (aggregator_real32_3d_t) this%aggregated_data = this%aggregated_data + ( & scale_val * this%source_data / div + offset_val - this%aggregated_data & ) / (this%counter + 1) type is (aggregator_real64_0d_t) this%aggregated_data = this%aggregated_data + ( & real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div + offset_val, kind=real32 & #else scale_val * this%source_data / div + offset_val, kind=real64 & #endif ) - this%aggregated_data & ) / (this%counter + 1) type is (aggregator_real64_1d_t) this%aggregated_data = this%aggregated_data + ( & real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div + offset_val, kind=real32 & #else scale_val * this%source_data / div + offset_val, kind=real64 & #endif ) - this%aggregated_data & ) / (this%counter + 1) type is (aggregator_real64_2d_t) this%aggregated_data = this%aggregated_data + ( & real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div + offset_val, kind=real32 & #else scale_val * this%source_data / div + offset_val, kind=real64 & #endif ) - this%aggregated_data & ) / (this%counter + 1) type is (aggregator_real64_3d_t) this%aggregated_data = this%aggregated_data + ( & real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div + offset_val, kind=real32 & #else scale_val * this%source_data / div + offset_val, kind=real64 & #endif ) - this%aggregated_data & ) / (this%counter + 1) class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select this%counter = this%counter + 1 end subroutine mean_accumulate subroutine sum_accumulate(this, scale, div, offset) !* Accumulate the aggregated data from the source data using the sum ! aggregation method. class(aggregator_t), intent(inout) :: this real, intent(in), optional :: scale !* An optional scaling factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: div !* An optional division factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: offset !* An optional offset to add to the source data before accumulation. ! Defaults to 0.0 if not provided. real :: scale_val, div_val, offset_val call get_accumulate_args(scale, div, offset, scale_val, div_val, offset_val) select type (this) type is (aggregator_int32_0d_t) this%aggregated_data = this%aggregated_data + int( & scale_val * this%source_data / div_val + offset_val & ) type is (aggregator_int32_1d_t) this%aggregated_data = this%aggregated_data + int( & scale_val * this%source_data / div_val + offset_val & ) type is (aggregator_int32_2d_t) this%aggregated_data = this%aggregated_data + int( & scale_val * this%source_data / div_val + offset_val & ) type is (aggregator_int32_3d_t) this%aggregated_data = this%aggregated_data + int( & scale_val * this%source_data / div_val + offset_val & ) type is (aggregator_real32_0d_t) this%aggregated_data = this%aggregated_data + ( & scale_val * this%source_data / div_val + offset_val & ) type is (aggregator_real32_1d_t) this%aggregated_data = this%aggregated_data + ( & scale_val * this%source_data / div_val + offset_val & ) type is (aggregator_real32_2d_t) this%aggregated_data = this%aggregated_data + ( & scale_val * this%source_data / div_val + offset_val & ) type is (aggregator_real32_3d_t) this%aggregated_data = this%aggregated_data + ( & scale_val * this%source_data / div_val + offset_val & ) type is (aggregator_real64_0d_t) this%aggregated_data = this%aggregated_data + real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif ) type is (aggregator_real64_1d_t) this%aggregated_data = this%aggregated_data + real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif ) type is (aggregator_real64_2d_t) this%aggregated_data = this%aggregated_data + real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif ) type is (aggregator_real64_3d_t) this%aggregated_data = this%aggregated_data + real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif ) class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select this%counter = this%counter + 1 end subroutine sum_accumulate subroutine point_accumulate(this, scale, div, offset) !* Accumulate the aggregated data from the source data using the point ! aggregation method. class(aggregator_t), intent(inout) :: this real, intent(in), optional :: scale !* An optional scaling factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: div !* An optional division factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: offset !* An optional offset to add to the source data before accumulation. ! Defaults to 0.0 if not provided. real :: scale_val, div_val, offset_val call get_accumulate_args(scale, div, offset, scale_val, div_val, offset_val) select type (this) type is (aggregator_int32_0d_t) this%aggregated_data = int(scale_val * this%source_data / div_val + offset_val) type is (aggregator_int32_1d_t) this%aggregated_data = int(scale_val * this%source_data / div_val + offset_val) type is (aggregator_int32_2d_t) this%aggregated_data = int(scale_val * this%source_data / div_val + offset_val) type is (aggregator_int32_3d_t) this%aggregated_data = int(scale_val * this%source_data / div_val + offset_val) type is (aggregator_real32_0d_t) this%aggregated_data = scale_val * this%source_data / div_val + offset_val type is (aggregator_real32_1d_t) this%aggregated_data = scale_val * this%source_data / div_val + offset_val type is (aggregator_real32_2d_t) this%aggregated_data = scale_val * this%source_data / div_val + offset_val type is (aggregator_real32_3d_t) this%aggregated_data = scale_val * this%source_data / div_val + offset_val type is (aggregator_real64_0d_t) this%aggregated_data = real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif ) type is (aggregator_real64_1d_t) this%aggregated_data = real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif ) type is (aggregator_real64_2d_t) this%aggregated_data = real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif ) type is (aggregator_real64_3d_t) this%aggregated_data = real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif ) class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select this%counter = this%counter + 1 end subroutine point_accumulate subroutine min_accumulate(this, scale, div, offset) !* Accumulate the aggregated data from the source data using the min ! aggregation method. class(aggregator_t), intent(inout) :: this real, intent(in), optional :: scale !* An optional scaling factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: div !* An optional division factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: offset !* An optional offset to add to the source data before accumulation. ! Defaults to 0.0 if not provided. real :: scale_val, div_val, offset_val call get_accumulate_args(scale, div, offset, scale_val, div_val, offset_val) select type (this) type is (aggregator_int32_0d_t) this%aggregated_data = min(this%aggregated_data, int( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_int32_1d_t) this%aggregated_data = min(this%aggregated_data, int( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_int32_2d_t) this%aggregated_data = min(this%aggregated_data, int( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_int32_3d_t) this%aggregated_data = min(this%aggregated_data, int( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_real32_0d_t) this%aggregated_data = min(this%aggregated_data, ( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_real32_1d_t) this%aggregated_data = min(this%aggregated_data, ( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_real32_2d_t) this%aggregated_data = min(this%aggregated_data, ( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_real32_3d_t) this%aggregated_data = min(this%aggregated_data, ( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_real64_0d_t) this%aggregated_data = min(this%aggregated_data, real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif )) type is (aggregator_real64_1d_t) this%aggregated_data = min(this%aggregated_data, real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif )) type is (aggregator_real64_2d_t) this%aggregated_data = min(this%aggregated_data, real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif )) type is (aggregator_real64_3d_t) this%aggregated_data = min(this%aggregated_data, real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif )) class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select this%counter = this%counter + 1 end subroutine min_accumulate subroutine max_accumulate(this, scale, div, offset) !* Accumulate the aggregated data from the source data using the max ! aggregation method. class(aggregator_t), intent(inout) :: this real, intent(in), optional :: scale !* An optional scaling factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: div !* An optional division factor to apply to the source data before ! accumulation. Defaults to 1.0 if not provided. real, intent(in), optional :: offset !* An optional offset to add to the source data before accumulation. ! Defaults to 0.0 if not provided. real :: scale_val, div_val, offset_val call get_accumulate_args(scale, div, offset, scale_val, div_val, offset_val) select type (this) type is (aggregator_int32_0d_t) this%aggregated_data = max(this%aggregated_data, int( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_int32_1d_t) this%aggregated_data = max(this%aggregated_data, int( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_int32_2d_t) this%aggregated_data = max(this%aggregated_data, int( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_int32_3d_t) this%aggregated_data = max(this%aggregated_data, int( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_real32_0d_t) this%aggregated_data = max(this%aggregated_data, ( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_real32_1d_t) this%aggregated_data = max(this%aggregated_data, ( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_real32_2d_t) this%aggregated_data = max(this%aggregated_data, ( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_real32_3d_t) this%aggregated_data = max(this%aggregated_data, ( & scale_val * this%source_data / div_val + offset_val & )) type is (aggregator_real64_0d_t) this%aggregated_data = max(this%aggregated_data, real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif )) type is (aggregator_real64_1d_t) this%aggregated_data = max(this%aggregated_data, real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif )) type is (aggregator_real64_2d_t) this%aggregated_data = max(this%aggregated_data, real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif )) type is (aggregator_real64_3d_t) this%aggregated_data = max(this%aggregated_data, real( & #ifdef ENFORCE_SINGLE_PRECISION scale_val * this%source_data / div_val + offset_val, kind=real32 & #else scale_val * this%source_data / div_val + offset_val, kind=real64 & #endif )) class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select this%counter = this%counter + 1 end subroutine max_accumulate subroutine point_reset(this) !* Reset the aggregated data for the point aggregation method. This is a ! no-op since point aggregation always takes the value of the most recent data ! point. class(aggregator_t), intent(inout) :: this end subroutine point_reset subroutine min_reset(this) !! Reset the aggregated data for the min aggregation method. class(aggregator_t), intent(inout) :: this select type (this) type is (aggregator_int32_0d_t) this%aggregated_data = huge(0_int32) type is (aggregator_int32_1d_t) this%aggregated_data = huge(0_int32) type is (aggregator_int32_2d_t) this%aggregated_data = huge(0_int32) type is (aggregator_int32_3d_t) this%aggregated_data = huge(0_int32) type is (aggregator_real32_0d_t) this%aggregated_data = huge(0.0_real32) type is (aggregator_real32_1d_t) this%aggregated_data = huge(0.0_real32) type is (aggregator_real32_2d_t) this%aggregated_data = huge(0.0_real32) type is (aggregator_real32_3d_t) this%aggregated_data = huge(0.0_real32) type is (aggregator_real64_0d_t) #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data = huge(0.0_real32) #else this%aggregated_data = huge(0.0_real64) #endif type is (aggregator_real64_1d_t) #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data = huge(0.0_real32) #else this%aggregated_data = huge(0.0_real64) #endif type is (aggregator_real64_2d_t) #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data = huge(0.0_real32) #else this%aggregated_data = huge(0.0_real64) #endif type is (aggregator_real64_3d_t) #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data = huge(0.0_real32) #else this%aggregated_data = huge(0.0_real64) #endif class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select this%counter = 0 end subroutine min_reset subroutine max_reset(this) !! Reset the aggregated data for the max aggregation method. class(aggregator_t), intent(inout) :: this select type (this) type is (aggregator_int32_0d_t) this%aggregated_data = -huge(0_int32) type is (aggregator_int32_1d_t) this%aggregated_data = -huge(0_int32) type is (aggregator_int32_2d_t) this%aggregated_data = -huge(0_int32) type is (aggregator_int32_3d_t) this%aggregated_data = -huge(0_int32) type is (aggregator_real32_0d_t) this%aggregated_data = -huge(0.0_real32) type is (aggregator_real32_1d_t) this%aggregated_data = -huge(0.0_real32) type is (aggregator_real32_2d_t) this%aggregated_data = -huge(0.0_real32) type is (aggregator_real32_3d_t) this%aggregated_data = -huge(0.0_real32) type is (aggregator_real64_0d_t) #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data = -huge(0.0_real32) #else this%aggregated_data = -huge(0.0_real64) #endif type is (aggregator_real64_1d_t) #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data = -huge(0.0_real32) #else this%aggregated_data = -huge(0.0_real64) #endif type is (aggregator_real64_2d_t) #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data = -huge(0.0_real32) #else this%aggregated_data = -huge(0.0_real64) #endif type is (aggregator_real64_3d_t) #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data = -huge(0.0_real32) #else this%aggregated_data = -huge(0.0_real64) #endif class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select this%counter = 0 end subroutine max_reset subroutine other_reset(this) !! Reset the aggregated data for aggregation methods other than point, min, and max. class(aggregator_t), intent(inout) :: this select type (this) type is (aggregator_int32_0d_t) this%aggregated_data = 0_int32 type is (aggregator_int32_1d_t) this%aggregated_data = 0_int32 type is (aggregator_int32_2d_t) this%aggregated_data = 0_int32 type is (aggregator_int32_3d_t) this%aggregated_data = 0_int32 type is (aggregator_real32_0d_t) this%aggregated_data = 0.0_real32 type is (aggregator_real32_1d_t) this%aggregated_data = 0.0_real32 type is (aggregator_real32_2d_t) this%aggregated_data = 0.0_real32 type is (aggregator_real32_3d_t) this%aggregated_data = 0.0_real32 type is (aggregator_real64_0d_t) #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data = 0.0_real32 #else this%aggregated_data = 0.0_real64 #endif type is (aggregator_real64_1d_t) #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data = 0.0_real32 #else this%aggregated_data = 0.0_real64 #endif type is (aggregator_real64_2d_t) #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data = 0.0_real32 #else this%aggregated_data = 0.0_real64 #endif type is (aggregator_real64_3d_t) #ifdef ENFORCE_SINGLE_PRECISION this%aggregated_data = 0.0_real32 #else this%aggregated_data = 0.0_real64 #endif class default call cable_abort("Unexpected aggregator type.", file=__FILE__, line=__LINE__) end select this%counter = 0 end subroutine other_reset function new_aggregator_int32_0d_t(source_data) result(agg) !! Create a new 0D integer aggregator. integer(kind=int32), intent(inout), target :: source_data !! The source data array to be sampled by the aggregator. type(aggregator_int32_0d_t) :: agg agg%source_data => source_data end function new_aggregator_int32_0d_t function new_aggregator_int32_1d_t(source_data) result(agg) !! Create a new 1D integer aggregator. integer(kind=int32), dimension(:), intent(inout), target :: source_data !! The source data array to be sampled by the aggregator. type(aggregator_int32_1d_t) :: agg agg%source_data => source_data end function new_aggregator_int32_1d_t function new_aggregator_int32_2d_t(source_data) result(agg) !! Create a new 2D integer aggregator. integer(kind=int32), dimension(:,:), intent(inout), target :: source_data !! The source data array to be sampled by the aggregator. type(aggregator_int32_2d_t) :: agg agg%source_data => source_data end function new_aggregator_int32_2d_t function new_aggregator_int32_3d_t(source_data) result(agg) !! Create a new 3D integer aggregator. integer(kind=int32), dimension(:,:,:), intent(inout), target :: source_data !! The source data array to be sampled by the aggregator. type(aggregator_int32_3d_t) :: agg agg%source_data => source_data end function new_aggregator_int32_3d_t function new_aggregator_real32_0d(source_data) result(agg) !! Create a new 0D 32-bit real aggregator. real(kind=real32), intent(inout), target :: source_data !! The source data array to be sampled by the aggregator. type(aggregator_real32_0d_t) :: agg agg%source_data => source_data end function new_aggregator_real32_0d function new_aggregator_real32_1d(source_data) result(agg) !! Create a new 1D 32-bit real aggregator. real(kind=real32), dimension(:), intent(inout), target :: source_data !! The source data array to be sampled by the aggregator. type(aggregator_real32_1d_t) :: agg agg%source_data => source_data end function new_aggregator_real32_1d function new_aggregator_real32_2d(source_data) result(agg) !! Create a new 2D 32-bit real aggregator. real(kind=real32), dimension(:,:), intent(inout), target :: source_data !! The source data array to be sampled by the aggregator. type(aggregator_real32_2d_t) :: agg agg%source_data => source_data end function new_aggregator_real32_2d function new_aggregator_real32_3d(source_data) result(agg) !! Create a new 3D 32-bit real aggregator. real(kind=real32), dimension(:,:,:), intent(inout), target :: source_data !! The source data array to be sampled by the aggregator. type(aggregator_real32_3d_t) :: agg agg%source_data => source_data end function new_aggregator_real32_3d function new_aggregator_real64_0d(source_data) result(agg) !! Create a new 0D 64-bit real aggregator. real(kind=real64), intent(inout), target :: source_data !! The source data array to be sampled by the aggregator. type(aggregator_real64_0d_t) :: agg agg%source_data => source_data end function new_aggregator_real64_0d function new_aggregator_real64_1d(source_data) result(agg) !! Create a new 1D 64-bit real aggregator. real(kind=real64), dimension(:), intent(inout), target :: source_data !! The source data array to be sampled by the aggregator. type(aggregator_real64_1d_t) :: agg agg%source_data => source_data end function new_aggregator_real64_1d function new_aggregator_real64_2d(source_data) result(agg) !! Create a new 2D 64-bit real aggregator. real(kind=real64), dimension(:,:), intent(inout), target :: source_data !! The source data array to be sampled by the aggregator. type(aggregator_real64_2d_t) :: agg agg%source_data => source_data end function new_aggregator_real64_2d function new_aggregator_real64_3d(source_data) result(agg) !! Create a new 3D 64-bit real aggregator. real(kind=real64), dimension(:,:,:), intent(inout), target :: source_data !! The source data array to be sampled by the aggregator. type(aggregator_real64_3d_t) :: agg agg%source_data => source_data end function new_aggregator_real64_3d end module