fix: improve the performance in nearest_index (NOAA-GFDL#1445)

axis_utils/include/axis_utils2.inc

@@ -306,64 +306,62 @@
     integer                               :: ia !< dimension of "array"
     integer                               :: i !< For looping through "array"
 
+    logical :: increasing !< .True. if the array is increasing
+
     ia = SIZE(array(:))
 
-    ! check if array is monotonous
+    ! check if array is increasing
+    increasing = .true.
     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
+       IF( array(i) .lt. array(i-1)) then
+         increasing = .false.
+         exit
+       endif
     END DO
 
-    if (array(1) < array(ia)) then
-    !< increasing array
+    if (.not. increasing) then
+      ! if not increasing, check that it is decreasing
+      DO i = 2, ia-1
+         IF( array(i) .gt. array(i-1)) &
+            call mpp_error(FATAL, 'axis_utils2::nearest_index array is NOT monotonously ordered')
+      END DO
+    endif
 
+    array_is_increasing: if (increasing) then
       !< Check if the rval is outside the range of the array
-      if (rval < array(1)) then
-         NEAREST_INDEX_ = 1
-         return
-      elseif (rval > array(ia)) then
-         NEAREST_INDEX_ = ia
-         return
+      if (rval .le. array(1)) then
+        NEAREST_INDEX_ = 1
+        return
+      elseif (rval .ge. array(ia)) then
+        NEAREST_INDEX_ = ia
+        return
       endif
 
-      DO i = 1, ia-1
-         IF ( (array(i)<=rval).AND.(array(i+1)>= rval) ) THEN
-            IF( rval - array(i) <= array(i+1) - rval ) THEN
-               NEAREST_INDEX_ = i
-               return
-            ELSE
-               NEAREST_INDEX_ = i+1
-               return
-            ENDIF
-            EXIT
-         END IF
+      DO i = 2, ia
+        if (rval .le. array(i)) then
+          NEAREST_INDEX_ = i
+          if (array(i) -rval .gt. rval - array(i-1)) NEAREST_INDEX_ = i - 1
+          return
+        endif
       END DO
-    else
-    !< Decreasing Array
-
+    else !array_is_decreasing
       !< Check if the rval is outside the range of the array
-      if (rval < array(ia)) then
-         NEAREST_INDEX_ = ia
-         return
-      elseif (rval > array(1)) then
-         NEAREST_INDEX_ = 1
-         return
+      if (rval .le. array(ia)) then
+        NEAREST_INDEX_ = ia
+        return
+      elseif (rval .gt. array(1)) then
+        NEAREST_INDEX_ = 1
+        return
       endif
 
-      DO i = 1, ia-1
-         IF ( (array(i)>=rval).AND.(array(i+1)<= rval) ) THEN
-            IF ( array(i)-rval <= rval-array(i+1) ) THEN
-               NEAREST_INDEX_ = i
-               return
-            ELSE
-               NEAREST_INDEX_ = i+1
-               return
-            END IF
-         END IF
+      DO i = 2, ia
+        if (rval .ge. array(i)) then
+          NEAREST_INDEX_ = i
+          if (rval - array(i) .gt. array(i-1) -rval ) NEAREST_INDEX_ = i - 1
+            return
+        endif
       END DO
-    endif
+    endif array_is_increasing
   end function NEAREST_INDEX_
 
   !#############################################################################

test_fms/axis_utils/test_axis_utils.F90

@@ -395,19 +395,19 @@ subroutine test_frac_index_fail
 
 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)
+  real(k) :: arr(7)
+  integer :: ans(16)
 
   if (increasing_array) then
-    arr = [5._k, 12._k, 20._k, 40._k, 100._k]
-    ans=(/1, 5, 1, 2, 3, 4, 5, 1, 2, 2, 3, 3/)
+    arr = [-6._k, -3._k, 5._k, 12._k, 20._k, 40._k, 100._k]
+    ans=(/1, 7, 3, 4, 5, 6, 7, 3, 4, 4, 5, 5, 1, 2, 1, 2/)
   else
-    arr = [100._k, 40._k, 20._k, 12._k, 5._k]
-    ans=(/5, 1, 5, 4, 3, 2, 1, 5, 4, 4, 3, 3/)
+    arr = [100._k, 40._k, 20._k, 12._k, 5._k, -3._k, -6._k]
+    ans=(/7, 1, 5, 4, 3, 2, 1, 5, 4, 4, 3, 3, 7, 6, 7, 6/)
   endif
 
   ! Test values beyond array boundaries
-  call nearest_index_assert(4._k,    arr, ans(1))
+  call nearest_index_assert(-7._k,    arr, ans(1))
   call nearest_index_assert(1000._k, arr, ans(2))
 
   ! Test values actually in the array
@@ -423,6 +423,13 @@ subroutine test_nearest_index(increasing_array)
   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))
+
+  ! Test the negative numbers
+  call nearest_index_assert(-6._k,    arr, ans(13))
+  call nearest_index_assert(-3._k,    arr, ans(14))
+  call nearest_index_assert(-5._k,    arr, ans(15))
+  call nearest_index_assert(-1._k,    arr, ans(16))
+
 end subroutine
 
 subroutine nearest_index_assert(val, arr, ret_expected)