update nearest_index to support motonically increasing/decreasing arrays

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).AND.array(i)>array(i+1)) .OR. (array(i-1)>array(i).AND.array(i)<array(i+1))) THEN
+          ! <ERROR STATUS="FATAL">array NOT monotonously ordered</ERROR>
+          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_
 

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)

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.