Updating interval division to use linear

Also improving the performance of linear interval division, when used for
large time intervals.

src/framework/mpas_timekeeping.F

@@ -1511,13 +1511,11 @@ subroutine mpas_interval_division(ref_time, num, den, n, rem)
       integer :: days, secondsNum, secondsDen
       integer (kind=I8KIND) :: seconds
 
-      call mpas_get_timeInterval(num, startTimeIn=ref_time, DD=days, S_i8=seconds, S_n=secondsNum, S_d=secondsDen)
-      call mpas_set_timeInterval(newNum, DD=days, S_i8=seconds, S_n=secondsNum, S_d=secondsDen)
-
-      call mpas_get_timeInterval(den, startTimeIn=ref_time, DD=days, S_i8=seconds, S_n=secondsNum, S_d=secondsDen)
-      call mpas_set_timeInterval(newDen, DD=days, S_i8=seconds, S_n=secondsNum, S_d=secondsDen)
-
-      call mpas_interval_division_log(newNum, newDen, n, rem)
+      if ( num % ti % YR == 0 .and. num % ti % MM == 0 .and. den % ti % YR == 0 .and. den % ti % MM == 0 ) then
+          call mpas_interval_division_log(num, den, n, rem)
+      else
+          call mpas_interval_division_linear(ref_time, num, den, n, rem)
+      end if
 
    end subroutine mpas_interval_division
 
@@ -1624,10 +1622,12 @@ subroutine mpas_interval_division_linear(ref_time, num, den, n, rem)
       integer, intent(out) :: n
       type (MPAS_TimeInterval_type), intent(out) :: rem
 
+      integer :: m
+
       type (MPAS_Time_type) :: target_time
-      type (MPAS_Time_type) :: updated_time
+      type (MPAS_Time_type) :: updated_time, mid_time
 
-      type (MPAS_TimeInterval_type) :: temp
+      type (MPAS_TimeInterval_type) :: temp, mid_int
       type (MPAS_TimeInterval_type) :: zero
 
       target_time = ref_time + num
@@ -1645,10 +1645,40 @@ subroutine mpas_interval_division_linear(ref_time, num, den, n, rem)
 
       ! One interval of den already fits into num
       n = n + 1
+      temp = den
+
+      ! Search forward, doubling the interval each time.
+      do while (target_time > updated_time)
+         n = n * 2
+         temp = den * n
+         updated_time = ref_time + temp
+      end do
+
+      ! Setup midpoint of search
+      ! The last value of n puts updated_time after target_time, need to back off and find the final time.
+      n = n / 2
+      m = n
+      mid_int = den * n
+      temp = mid_int
+      updated_time = ref_time + mid_int + temp
+
+      ! Seach backward, halving the interval each time.
+      do while (target_time < updated_time)
+         m  = m / 2
+         temp = den * m
+         updated_time = ref_time + mid_int + temp
+      end do
+
+      ! Final number of interavls is n + m
+      n = n + m
+
+      ! Do a final linear search, just to ensure we aren't missing any divisions.
+      temp = den * n
+      updated_time = ref_time + temp
 
       do while (target_time > updated_time)
-         updated_time = updated_time + den
          n = n + 1
+         updated_time = updated_time + den
       end do
 
       ! Here, if updated_time is larger than target time. Need to subtract den once, and compute remainder