feat: support monotonically increasing/decreasing array in data_overr…

…ide (#1388)

axis_utils/include/axis_utils2.inc

@@ -46,6 +46,9 @@
                                               !! 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) !< Index to use when reading the edges from the file
+                           !! (/1, 2/) if the axis data is monotonically increasing
+                           !! (/2, 1/) if the axis data is monotonically decreasing
 
   ndims = get_variable_num_dimensions(fileobj, name)
   allocate(dim_sizes(ndims))
@@ -108,8 +111,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)
@@ -265,7 +270,7 @@
   !!     inputs:
   !!
   !!     rval   = arbitrary data...same units as elements in "array"
-  !!     array  = array of data points  (must be monotonically increasing)
+  !!     array  = array of data points  (must be monotonically)
   !!     ia     = dimension of "array"
   !!
   !!     output:
@@ -293,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_(rval, array)
+    real(kind=FMS_AU_KIND_), intent(in)               :: rval !< arbitrary data...same units as elements in "array"
+    real(kind=FMS_AU_KIND_), intent(in), dimension(:) :: array !< array of data points  (must be monotonic)
 
     integer                               :: NEAREST_INDEX_
     integer                               :: ia !< dimension of "array"
-    integer                               :: i, ii, iunit
-    real(kind=FMS_AU_KIND_)               :: rval !< 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
-          iunit = stdout()
-          write (iunit,*) '=> Error: "nearest_index" array must be monotonically increasing' &
-                         & // 'when searching for nearest value to ', rval
-          write (iunit,*) '          array(i) < array(i-1) for i=',i
-          write (iunit,*) '          array(i) for i=1..ia follows:'
-          do ii = 1, ia
-             write (iunit,*) 'i=',ii, ' array(i)=', array(ii)
-          enddo
-          call mpp_error(FATAL,' "nearest_index" array must be monotonically increasing.')
-       endif
-    enddo
+    integer                               :: i !< For looping through "array"
+
+    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 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
+      endif
 
-    if (rval < array(1) .or. rval > array(ia)) then
-       if (rval < array(1))  NEAREST_INDEX_ = 1
-       if (rval > array(ia)) NEAREST_INDEX_ = ia
+      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
+      END DO
     else
-       i = 1
-       keep_going = .true.
-       do while (i <= ia .and. keep_going)
-          i = i+1
-          if (rval <= array(i)) then
-             NEAREST_INDEX_ = i
-             if (array(i)-rval > rval-array(i-1)) NEAREST_INDEX_ = i-1
-             keep_going = .false.
-          endif
-       enddo
+    !< Decreasing Array
+
+      !< 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
+      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
+      END DO
     endif
   end function NEAREST_INDEX_
 

data_override/include/data_override.inc

@@ -1279,14 +1279,14 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d
         ie_src = axis_sizes(1)
         js_src = 1
         je_src = axis_sizes(2)
-        do j = 1, axis_sizes(2)+1
-           if( override_array(curr_position)%lat_in(j) > lat_min ) exit
-           js_src = j
-        enddo
-        do j = 1, axis_sizes(2)+1
-           je_src = j
-           if( override_array(curr_position)%lat_in(j) > lat_max ) exit
-        enddo
+        ! do j = 1, axis_sizes(2)+1
+        !    if( override_array(curr_position)%lat_in(j) > lat_min ) exit
+        !    js_src = j
+        ! enddo
+        ! do j = 1, axis_sizes(2)+1
+        !    je_src = j
+        !    if( override_array(curr_position)%lat_in(j) > lat_max ) exit
+        ! enddo
 
         !--- bicubic interpolation need one extra point in each direction. Also add
         !--- one more point for because lat_in is in the corner but the interpolation

horiz_interp/horiz_interp_bilinear.F90

@@ -34,6 +34,7 @@ module horiz_interp_bilinear_mod
   use constants_mod,         only: PI
   use horiz_interp_type_mod, only: horiz_interp_type, stats
   use platform_mod,          only: r4_kind, r8_kind
+  use axis_utils2_mod,       only: nearest_index
 
   implicit none
   private
@@ -63,12 +64,6 @@ module horiz_interp_bilinear_mod
   integer, parameter :: DUMMY = -999
 
 !! Private helper routines, interfaces for mixed real precision support
-
-  interface indp
-    module procedure indp_r4
-    module procedure indp_r8
-  end interface
-
   interface intersect
     module procedure intersect_r4
     module procedure intersect_r8

horiz_interp/include/horiz_interp_bilinear.inc

@@ -34,12 +34,20 @@
     real(FMS_HI_KIND_)    :: wtw, wte, wts, wtn, lon, lat, tpi, hpi
     real(FMS_HI_KIND_)    :: glt_min, glt_max, gln_min, gln_max, min_lon, max_lon
     integer,parameter     :: kindl = FMS_HI_KIND_
+    logical               :: decreasing_lat !< .True. if latitude is monotically decreasing
+    logical               :: decreasing_lon !< .True. if longitude is monotically decreasing
 
     warns = 0
     if(present(verbose)) warns = verbose
     src_is_modulo = .true.
     if (present(src_modulo)) src_is_modulo = src_modulo
 
+    decreasing_lat = .false.
+    if (lat_in(1) > lat_in(2)) decreasing_lat = .true.
+
+    decreasing_lon = .false.
+    if (lon_in(1) > lon_in(2)) decreasing_lon = .true.
+
     hpi = 0.5_kindl * real(pi, FMS_HI_KIND_)
     tpi = 4.0_kindl * hpi
     glt_min = hpi
@@ -95,11 +103,23 @@
           glt_min = min(lat,glt_min);  glt_max = max(lat,glt_max)
           gln_min = min(lon,gln_min);  gln_max = max(lon,gln_max)
 
-          is = indp(lon, lon_in )
-          if( lon_in(is) .gt. lon ) is = max(is-1,1)
+          is = nearest_index(lon, lon_in )
+          if (decreasing_lon) then
+            ! Lon_in is increasing
+            ! This is so that is is the lower bound.
+            ! For example, if the array is [50 40 30 20 10] and lon is 11, `is` is going to be 5 from `nearest_index`
+            ! but it needs to be 4 so that it can use the data at lon_in(4) and lon_in(5)
+            if( lon_in(is) .lt. lon ) is = max(is-1,1)
+          else
+            ! Lon_in is increasing
+            ! This is so that is is the lower bound.
+            ! For example, if the array is [10 20 30 40 50] and lon is 49, `is` is going to be 5 from `nearest_index`
+            ! but it needs to be 4 so that it can use the data at lon_in(4) and lon_in(5)
+            if( lon_in(is) .gt. lon ) is = max(is-1,1)
+          endif
           if( lon_in(is) .eq. lon .and. is .eq. nlon_in) is = max(is - 1,1)
           ie = min(is+1,nlon_in)
-          if(lon_in(is) .ne. lon_in(ie) .and. lon_in(is) .le. lon) then
+          if(lon_in(is) .ne. lon_in(ie) .and. (decreasing_lon .or. lon_in(is) .le. lon)) then
              wtw = ( lon_in(ie) - lon) / (lon_in(ie) - lon_in(is) )
           else
              !     east or west of the last data value. this could be because a
@@ -115,13 +135,18 @@
           endif
           wte = 1.0_kindl - wtw
 
-          js = indp(lat, lat_in )
-
-          if( lat_in(js) .gt. lat ) js = max(js - 1, 1)
+          js = nearest_index(lat, lat_in )
+          if (decreasing_lat) then
+            ! Lat_in is decreasing
+            if( lat_in(js) .lt. lat ) js = max(js - 1, 1)
+          else
+            ! Lat_in is increasing
+            if( lat_in(js) .gt. lat ) js = max(js - 1, 1)
+          endif
           if( lat_in(js) .eq. lat .and. js .eq. nlat_in) js = max(js - 1, 1)
           je = min(js + 1, nlat_in)
 
-          if ( lat_in(js) .ne. lat_in(je) .and. lat_in(js) .le. lat) then
+          if ( lat_in(js) .ne. lat_in(je) .and. (decreasing_lat .or. lat_in(js) .le. lat)) then
              wts = ( lat_in(je) - lat )/(lat_in(je)-lat_in(js))
           else
              !     north or south of the last data value. this could be because a
@@ -1180,51 +1205,4 @@
     return
 
   end subroutine
-
-
-  !#######################################################################
-  !> @returns index of nearest data point to "value"
-  !! if "value" is outside the domain of "array" then INDP_ = 1
-  !! or "ia" depending on whether array(1) or array(ia) is
-  !! closest to "value"
-  function INDP_ (rval, array)
-    integer                        :: INDP_ !< index of nearest data point within "array"
-                                           !! corresponding to "value".
-    real(FMS_HI_KIND_), dimension(:), intent(in) :: array !< array of data points (must be monotonically increasing)
-    real(FMS_HI_KIND_), intent(in)               :: rval !< arbitrary data, same units as elements in 'array'
-
-    !=======================================================================
-
-    integer i, ia, iunit
-    logical keep_going
-    !
-    ia = size(array(:))
-    do i=2,ia
-       if (array(i) .lt. array(i-1)) then
-          iunit = stdout()
-          write (iunit,*) &
-               ' => Error: array must be monotonically increasing in "INDP_"' , &
-               '           when searching for nearest element to value=',rval
-          write (iunit,*) '           array(i) < array(i-1) for i=',i
-          write (iunit,*) '           array(i) for i=1..ia follows:'
-          call mpp_error()
-       endif
-    enddo
-    if (rval .lt. array(1) .or. rval .gt. array(ia)) then
-       if (rval .lt. array(1))  INDP_ = 1
-       if (rval .gt. array(ia)) INDP_ = ia
-    else
-       i=1
-       keep_going = .true.
-       do while (i .le. ia .and. keep_going)
-          i = i+1
-          if (rval .le. array(i)) then
-             INDP_ = i
-             if (array(i)-rval .gt. rval-array(i-1)) INDP_ = i-1
-             keep_going = .false.
-          endif
-       enddo
-    endif
-    return
-  end function INDP_
 !> @}

horiz_interp/include/horiz_interp_bilinear_r4.fh

@@ -45,8 +45,5 @@
 #undef INTERSECT_
 #define INTERSECT_ intersect_r4
 
-#undef INDP_
-#define INDP_ indp_r4
-
 #include "horiz_interp_bilinear.inc"
 !> @}

horiz_interp/include/horiz_interp_bilinear_r8.fh

@@ -45,8 +45,5 @@
 #undef INTERSECT_
 #define INTERSECT_ intersect_r8
 
-#undef INDP_
-#define INDP_ indp_r8
-
 #include "horiz_interp_bilinear.inc"
 !> @}