Support monotonically increasing/decreasing array in data_override

axis_utils/include/axis_utils2.inc

@@ -46,6 +46,7 @@
                                               !! indicating to reproduce the mpp_io bug where
                                               !! the null characters were not removed after reading a string attribute
   integer, parameter                                   :: lkind = FMS_AU_KIND_
+  integer :: edge_index(2)
 
   ndims = get_variable_num_dimensions(fileobj, name)
   allocate(dim_sizes(ndims))
@@ -108,8 +109,10 @@
 
       allocate(r2d(dim_sizes(1), dim_sizes(2)))
       call read_data(fileobj, buffer, r2d)
-      edge_data(1:dim_sizes(2)) = r2d(1,:)
-      edge_data(dim_sizes(2)+1) = r2d(2,dim_sizes(2))
+      edge_index = (/1, 2/)
+      if (r2d(1,1) .gt. r2d(1,2)) edge_index = (/2, 1 /)
+      edge_data(1:dim_sizes(2)) = r2d(edge_index(1),:)
+      edge_data(dim_sizes(2)+1) = r2d(edge_index(2),dim_sizes(2))
       deallocate(r2d)
     endif
     deallocate(dim_sizes)
@@ -267,7 +270,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 +298,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,17 +76,25 @@ 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)"
       call test_nearest_index_fail
 
-    case ('--axis-edges')
-      print "(A)", "Testing axis_edges"
-      call test_axis_edges
+    case ('--axis-edges-increasing')
+      print "(A)", "Testing axis_edges-increasing"
+      call test_axis_edges(.true.)
+
+    case ('--axis-edges-decreasing')
+      print "(A)", "Testing axis_edges-decreasing"
+      call test_axis_edges(.false.)
 
     case ('--tranlon')
       print "(A)", "Testing tranlon"
@@ -385,28 +393,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
+    arr = [5._k, 12._k, 20._k, 40._k, 100._k]
+    ans=(/1, 5, 1, 2, 3, 4, 5, 1, 2, 2, 3, 3/)
+  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/)
+  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)
@@ -433,24 +449,42 @@ subroutine test_nearest_index_fail
   ret_test = nearest_index(5._k, arr)
 end subroutine
 
-subroutine test_axis_edges
+subroutine test_axis_edges(increasing_array)
+  logical, intent(in) :: increasing_array !< .True. if test using an increasing array
   real(k) :: data_in_var(10)
   real(k) :: data_in_var_edges(2,10)
   real(k) :: data_in_answers(11)
   type(FmsNetcdfFile_t) :: fileobj
   real(k)    :: answers(11)
-  integer :: i
+  integer :: count
+  integer :: count_factor
+  integer :: factor
+
+  if (increasing_array) then
+    count = 0
+    factor = 1
+    count_factor = -1
+  else
+    count = 11
+    factor = -1
+    count_factor = 0
+  endif
 
   do i=1,10
-     data_in_var(i) = real(i, k) - 0.5_k
+     count = count + factor
+     data_in_var(i) = real(count, k) - 0.5_k
 
-     data_in_var_edges(1,i) = real(i-1, k)
-     data_in_var_edges(2,i) = real(i, k)
+     data_in_var_edges(1,i) = real(count-1, k)
+     data_in_var_edges(2,i) = real(count, k)
 
-     data_in_answers(i) = real(i-1, k)
+     data_in_answers(i) = real(count + count_factor, k)
   enddo
 
-  data_in_answers(11) = 10._k
+  if (increasing_array) then
+    data_in_answers(11) = real(count, k)
+  else
+    data_in_answers(11) = real(count + factor, k)
+  endif
 
   call open_netcdf_w(fileobj)
 

test_fms/axis_utils/test_axis_utils2.sh

@@ -27,11 +27,11 @@
 # 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-increasing --axis_edges-decreasing --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.
-SKIP_TESTS="test_axis_utils2.15 test_axis_utils2.16"
+SKIP_TESTS="test_axis_utils2.19 test_axis_utils2.20"
 
 # Run the tests