diff --git a/axis_utils/include/axis_utils2.inc b/axis_utils/include/axis_utils2.inc index baad2092e5..3943fd9c0a 100644 --- a/axis_utils/include/axis_utils2.inc +++ b/axis_utils/include/axis_utils2.inc @@ -267,7 +267,7 @@ !! inputs: !! !! value = arbitrary data...same units as elements in "array" - !! array = array of data points (must be monotonically increasing) + !! array = array of data points (must be monotonic) !! ia = dimension of "array" !! !! output: @@ -295,48 +295,71 @@ !! z(k1) would be the nearest data point to 12.5 and z(k2) would !! be the nearest data point to 0.0 !! @return integer nearest_index - - - function NEAREST_INDEX_(value, array) + real(kind=FMS_AU_KIND_) , intent(in) :: value !< arbitrary data...same units as elements in "array" + real(kind=FMS_AU_KIND_), dimension(:), intent(in) :: array !< array of data points (must be monotonic) integer :: NEAREST_INDEX_ integer :: ia !< dimension of "array" - integer :: i, ii, unit - real(kind=FMS_AU_KIND_) :: value !< arbitrary data...same units as elements in "array" - real(kind=FMS_AU_KIND_), dimension(:) :: array !< array of data points (must be monotonically increasing) - logical :: keep_going - - ia = size(array(:)) - - do i = 2, ia - if (array(i) < array(i-1)) then - unit = stdout() - write (unit,*) '=> Error: "nearest_index" array must be monotonically increasing' & - & // 'when searching for nearest value to ', value - write (unit,*) ' array(i) < array(i-1) for i=',i - write (unit,*) ' array(i) for i=1..ia follows:' - do ii = 1, ia - write (unit,*) 'i=',ii, ' array(i)=', array(ii) - enddo - call mpp_error(FATAL,' "nearest_index" array must be monotonically increasing.') - endif - enddo + integer :: i !< For do loops + + ia = SIZE(array(:)) + + ! check if array is monotonous + DO i = 2, ia-1 + IF( (array(i-1)array(i+1)) .OR. (array(i-1)>array(i).AND.array(i)array NOT monotonously ordered + CALL mpp_error(FATAL, 'axis_utils2::nearest_index array is NOT monotonously ordered') + END IF + END DO + + if (array(1) < array(ia)) then + !< increasing array + + !< Check if the value is outside the range of the array + if (value < array(1)) then + NEAREST_INDEX_ = 1 + return + elseif (value > array(ia)) then + NEAREST_INDEX_ = ia + return + endif - if (value < array(1) .or. value > array(ia)) then - if (value < array(1)) NEAREST_INDEX_ = 1 - if (value > array(ia)) NEAREST_INDEX_ = ia + DO i = 1, ia-1 + IF ( (array(i)<=value).AND.(array(i+1)>= value) ) THEN + IF( value - array(i) <= array(i+1) - value ) THEN + NEAREST_INDEX_ = i + return + ELSE + NEAREST_INDEX_ = i+1 + return + ENDIF + EXIT + END IF + END DO else - i = 1 - keep_going = .true. - do while (i <= ia .and. keep_going) - i = i+1 - if (value <= array(i)) then - NEAREST_INDEX_ = i - if (array(i)-value > value-array(i-1)) NEAREST_INDEX_ = i-1 - keep_going = .false. - endif - enddo + !< Decreasing Array + + !< Check if the value is outside the range of the array + if (value < array(ia)) then + NEAREST_INDEX_ = ia + return + elseif (value > array(1)) then + NEAREST_INDEX_ = 1 + return + endif + + DO i = 1, ia-1 + IF ( (array(i)>=value).AND.(array(i+1)<= value) ) THEN + IF ( array(i)-value <= value-array(i+1) ) THEN + NEAREST_INDEX_ = i + return + ELSE + NEAREST_INDEX_ = i+1 + return + END IF + END IF + END DO endif end function NEAREST_INDEX_ diff --git a/test_fms/axis_utils/test_axis_utils.F90 b/test_fms/axis_utils/test_axis_utils.F90 index aac74de010..c2765bae4b 100644 --- a/test_fms/axis_utils/test_axis_utils.F90 +++ b/test_fms/axis_utils/test_axis_utils.F90 @@ -76,9 +76,13 @@ program test_axis_utils print "(A)", "Testing frac_index (FAILURE)" call test_frac_index_fail - case ('--nearest-index') - print "(A)", "Testing nearest_index" - call test_nearest_index + case ('--nearest-index-increasing') + print "(A)", "Testing nearest_index with a monotonically increasing array" + call test_nearest_index(.true.) + + case ('--nearest-index-decreasing') + print "(A)", "Testing nearest_index with a monotonically decreasing array" + call test_nearest_index(.false.) case ('--nearest-index-fail') print "(A)", "Testing nearest_index (FAILURE)" @@ -385,28 +389,36 @@ subroutine test_frac_index_fail ret_test = frac_index(1.5_k, values) end subroutine -subroutine test_nearest_index +subroutine test_nearest_index(increasing_array) + logical, intent(in) :: increasing_array !< .True. if test using an increasing array real(k) :: arr(5) + integer :: ans(12) arr = [5._k, 12._k, 20._k, 40._k, 100._k] + if (increasing_array) then + ans=(/1, 5, 1, 2, 3, 4, 5, 1, 2, 2, 3, 3/) + else + arr = arr(ubound(arr,dim=1)::-1) + ans=(/5, 1, 5, 4, 3, 2, 1, 5, 4, 4, 3, 3/) + endif ! Test values beyond array boundaries - call nearest_index_assert(4._k, arr, 1) - call nearest_index_assert(1000._k, arr, size(arr)) + call nearest_index_assert(4._k, arr, ans(1)) + call nearest_index_assert(1000._k, arr, ans(2)) ! Test values actually in the array - call nearest_index_assert(5._k, arr, 1) - call nearest_index_assert(12._k, arr, 2) - call nearest_index_assert(20._k, arr, 3) - call nearest_index_assert(40._k, arr, 4) - call nearest_index_assert(100._k, arr, 5) + call nearest_index_assert(5._k, arr, ans(3)) + call nearest_index_assert(12._k, arr, ans(4)) + call nearest_index_assert(20._k, arr, ans(5)) + call nearest_index_assert(40._k, arr, ans(6)) + call nearest_index_assert(100._k, arr, ans(7)) ! Test the intervals between array values - call nearest_index_assert(6._k, arr, 1) - call nearest_index_assert(11._k, arr, 2) - call nearest_index_assert(15._k, arr, 2) - call nearest_index_assert(18._k, arr, 3) - call nearest_index_assert(29._k, arr, 3) + call nearest_index_assert(6._k, arr, ans(8)) + call nearest_index_assert(11._k, arr, ans(9)) + call nearest_index_assert(15._k, arr, ans(10)) + call nearest_index_assert(18._k, arr, ans(11)) + call nearest_index_assert(29._k, arr, ans(12)) end subroutine subroutine nearest_index_assert(val, arr, ret_expected) diff --git a/test_fms/axis_utils/test_axis_utils2.sh b/test_fms/axis_utils/test_axis_utils2.sh index 1288822481..eacfb89473 100755 --- a/test_fms/axis_utils/test_axis_utils2.sh +++ b/test_fms/axis_utils/test_axis_utils2.sh @@ -27,7 +27,7 @@ # Prepare the directory to run the tests. touch input.nml -TESTS_SUCCESS='--get-axis-modulo --get-axis-modulo-times --get-axis-cart --lon-in-range --frac-index --nearest-index --axis-edges --tranlon --interp-1d-1d --interp-1d-2d --interp-1d-3d' +TESTS_SUCCESS='--get-axis-modulo --get-axis-modulo-times --get-axis-cart --lon-in-range --frac-index --nearest-index-increasing --nearest-index-decreasing --axis-edges --tranlon --interp-1d-1d --interp-1d-2d --interp-1d-3d' TESTS_FAIL='--frac-index-fail --nearest-index-fail' # TODO: Enable these tests after tranlon's memory corruption bug is fixed.