!-----------------------------------------------------------------------------
!
! find_runs.f90 -- Generic routine to find runs of almost-same values
!		   in 1-D data series.
!
! 2017-apr-11	Original version for bias correction.
!		By Dave Allured, NOAA/ESRL/PSD/CIRES.
!
! This version finds multiple runs of almost-same values.
! Missing values are not included, and will interrupt runs.
!
! The result is a logical mask of the same length as the input
! array, and True in all positions of all qualifying runs.
!
! It is possible to have overlapping runs with slightly differing
! match tolerances.  Such runs will appear as a single long run
! in the result mask.
!
!-----------------------------------------------------------------------------

module find__runs
contains

subroutine find_runs (x, vmiss, min_length, tolerance_rel, mask_out)
   use config, only : dp
   implicit none

   real(dp), intent (in) :: x(:)	     ! 1-D input data series
   real(dp), intent (in) :: vmiss	     ! missing value in data series
   integer,  intent (in) :: min_length	     ! minimum allowable run length
   real(dp), intent (in) :: tolerance_rel    ! relative tolerance to match data

   logical,  intent(out), allocatable :: mask_out(:)  ! result, true within runs

! Local variables.

   integer i, j, k, ndata
   real(dp) thresh_abs

! Initialize.

   ndata = size (x)

   allocate (mask_out(ndata))
   mask_out(:) = .false.		! init result mask, false = no runs

! Use the three-loop method (i, j, k) because it is simpler to code
! and understand.

! Loop to check EVERY data point as the possible start of a run,
! up to the last possible starting point.

iloop: &
   do i = 1, ndata - min_length + 1
      if (x(i) == vmiss) cycle iloop	! missing value disqualifies a run

! Search for a run of the specified minimum length, within the given tolerance.
! Use an absolute threshold, not relative, to protect from divide by zero.

      thresh_abs = abs (tolerance_rel * x(i))

jloop: &
      do j = i + 1, i + min_length - 1
         if (x(j) == vmiss) cycle iloop    ! missing value disqualifies a run
         if (abs (x(j) - x(i)) > thresh_abs) cycle iloop   !  disqualify if
      end do jloop			       ! mismatch within minimum length

! Found a run of minimum qualifying length.  Now find the end of the run.

kloop: &
      do k = i + min_length - 1, ndata	  ! first iteration guaranteed to match
         if (x(k) == vmiss) exit kloop	  ! missing value terminates a run
         if (abs (x(k) - x(i)) > thresh_abs) exit kloop   ! mismatch terminates
         j = k				  ! j saves the highest matching index
      end do kloop

! Now found a complete run of near-constant values.
! Indicate all values within run.

      mask_out(i:j) = .true.		! true = values within current run
   end do iloop

end subroutine find_runs
end module find__runs