module array_utils ! ------------------------------------------------------------------------ ! !DESCRIPTION: ! This module contains routines for working with arrays. ! ! NOTES ! Subroutines transpose_wrapper_* and pack_wrapper could be moved to csm_share ! if they stick around for long; they are just here to work around the fact ! that we don't specify -Mallocatable=03 with pgi (which in turn is because ! that flag led to bugs with pgi14). ! ! !USES: #include "shr_assert.h" use shr_kind_mod, only : r8 => shr_kind_r8, r4 => shr_kind_r4, i4=>shr_kind_i4 use shr_log_mod , only : errMsg => shr_log_errMsg use shr_sys_mod , only : shr_sys_abort implicit none private save ! Public routines public :: find_k_max_indices ! returns the indices of the maximum k values in data, in sorted order public :: convert_to_logical ! convert a numeric array to logical array public :: transpose_wrapper ! wrap the intrinsic transpose function, first allocating the destination array public :: pack_wrapper ! wrap the intrinsic pack function, first allocating the destination array interface convert_to_logical !DIMS 1,2,3 !TYPE double,real,int module procedure convert_to_logical_{DIMS}d_{TYPE} end interface convert_to_logical interface transpose_wrapper !TYPE int,double module procedure transpose_wrapper_{TYPE} end interface transpose_wrapper character(len=*), parameter, private :: sourcefile = & __FILE__ contains !----------------------------------------------------------------------- subroutine find_k_max_indices(data, lb, k, max_indices) ! ! !DESCRIPTION: ! Returns the indices of the maximum k values in data, in sorted order (so the first ! index will be the largest, the second index will be the second largest, etc.) ! ! The output array (max_indices) should be of size k ! ! Speed estimate: O(n*k) (where n is the size of the input array) ! ! !ARGUMENTS: integer , intent(in) :: lb ! lower bound of data array real(r8) , intent(in) :: data(lb:) integer , intent(in) :: k ! number of max indices to find integer , intent(out) :: max_indices(:) ! ! !LOCAL VARIABLES: real(r8) :: max_vals(k) integer :: i, j integer :: insertion_location character(len=*), parameter :: subname = 'find_k_max_indices' !----------------------------------------------------------------------- SHR_ASSERT_FL((size(max_indices) == k), subname, __LINE__) if (k < 1 .or. k > size(data)) then call shr_sys_abort(subname//': must have 1 <= k <= size(data)') end if max_indices(:) = lb - 1 max_vals(:) = -huge(1._r8) do i = lbound(data,1), ubound(data,1) ! Determine where (if anywhere) data(i) should go in max_vals, relative to the ! values we've seen so far. Note that max_vals will be arranged in sorted order, ! with the largest value first. insertion_location = -1 do j = 1, k if (data(i) > max_vals(j)) then insertion_location = j exit end if end do if (insertion_location > 0) then ! Shift lower values right by one (note that the k'th value will be dropped) do j = k-1, insertion_location, -1 max_indices(j+1) = max_indices(j) max_vals(j+1) = max_vals(j) end do max_indices(insertion_location) = i max_vals(insertion_location) = data(i) end if end do end subroutine find_k_max_indices !----------------------------------------------------------------------- !DIMS 1,2,3 !TYPE double,real,int subroutine convert_to_logical_{DIMS}d_{TYPE}(data_numeric, data_logical) ! ! !DESCRIPTION: ! Convert a numeric variable to logical ! ! !ARGUMENTS: {VTYPE}, intent(in) :: data_numeric{DIMSTR} logical, intent(out) :: data_logical{DIMSTR} ! ! !LOCAL VARIABLES: character(len=*), parameter :: subname = 'convert_to_logical_{DIMS}d_{TYPE}' !----------------------------------------------------------------------- if (any(data_numeric /= 0 .and. data_numeric /= 1)) then call shr_sys_abort('convert_to_logical: bad value(s) for logical data') end if where (data_numeric == 1) data_logical = .true. elsewhere data_logical = .false. end where end subroutine convert_to_logical_{DIMS}d_{TYPE} !----------------------------------------------------------------------- !TYPE int,double subroutine transpose_wrapper_{TYPE}(arr_out, arr_in) ! ! !DESCRIPTION: ! Wrap the intrinsic transpose function, doing the necessary allocation of the ! destination array ! ! NOTE(wjs, 2015-10-21) This is supposed to be handled for you in Fortran2003. Within ! CESM, it currently appears to be handled for all compilers except pgi: It is ! handled with intel through the use of '-assume realloc_lhs'. It could be handled ! with pgi through the use of '-Mallocatable=03', but we currently don't use that ! flag, because it triggered bugs with pgi14. ! ! !USES: ! ! !ARGUMENTS: {VTYPE}, allocatable, intent(out) :: arr_out(:,:) {VTYPE}, intent(in) :: arr_in(:,:) ! ! !LOCAL VARIABLES: integer :: size1, size2 character(len=*), parameter :: subname = 'transpose_wrapper_{TYPE}' !----------------------------------------------------------------------- size1 = size(arr_in, 1) size2 = size(arr_in, 2) allocate(arr_out(size2, size1)) arr_out = transpose(arr_in) end subroutine transpose_wrapper_{TYPE} !----------------------------------------------------------------------- subroutine pack_wrapper(arr_packed, arr, mask) ! ! !DESCRIPTION: ! Wrap the intrinsic pack function, doing the necessary allocation of the destination ! array. ! ! NOTE(wjs, 2015-10-21) This is supposed to be handled for you in Fortran2003. Within ! CESM, it currently appears to be handled for all compilers except pgi: It is ! handled with intel through the use of '-assume realloc_lhs'. It could be handled ! with pgi through the use of '-Mallocatable=03', but we currently don't use that ! flag, because it triggered bugs with pgi14. ! ! !USES: ! ! !ARGUMENTS: real(r8), allocatable, intent(out) :: arr_packed(:) real(r8), intent(in) :: arr(:) logical , intent(in) :: mask(:) ! ! !LOCAL VARIABLES: integer :: packed_size character(len=*), parameter :: subname = 'pack_wrapper' !----------------------------------------------------------------------- packed_size = count(mask) allocate(arr_packed(packed_size)) arr_packed = pack(arr, mask) end subroutine pack_wrapper end module array_utils