! $Id: remove_spikes.f90 2022-10-26 09:51:28Z ychen $

!****s* fsi/remove_spikes *
!
! NAME
!   remove_spikes  - remove outliers in a timeseries
!
! SYNOPSIS
!    call remove_spikes(xin, yin, nx, threshold, xout, yout, idx, nout)
!
! DESCRIPTION
!    This subroutine to remove outliers in a timeseries based on segment mean and
!                    standard deviation
! INPUTS
!    REAL(wp), DIMENSION(:),    :: xin          absissa value (only real values)
!    REAL(wp), DIMENSION(:),    :: yin          ordinate value (only real values)
!    INTEGER,                   :: nx           segment length (nx << nmax and nx must be odd integer)
!    REAL(wp),                  :: threshold    threshold value to determine cutoff
!
! OUTPUT
!    REAL(wp), DIMENSION(:),    :: xout         absissa value (size = nout)
!    REAL(wp), DIMENSION(:),    :: yout         ordinate value (size = nout)
!    INTEGER,  DIMENSION(:),    :: idx          index of used values
!    INTEGER,                   :: nout         size of output data
!
! AUTHOR
!   Yong Chen, NOAA/NESDIS/STAR, yong.chen@noaa.gov
!   Loknath Adhikari, UMD/CISESS, loknath.adhikari@noaa.gov
!
! COPYRIGHT
!   Copyright (c) 2022-2023 Yong Chen <yong.chen@noaa.gov>
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

SUBROUTINE remove_spikes(xin, yin, nx, threshold, xout, yout, idx, nout)

      USE typesizes,         ONLY: wp => EightByteReal

      IMPLICIT none
      REAL(wp), DIMENSION(:),   INTENT(in)    :: xin
      REAL(wp), DIMENSION(:),   INTENT(in)    :: yin
      INTEGER,                  INTENT(inout) :: nx
      REAL(wp),                 INTENT(in)    :: threshold
      REAL(wp), DIMENSION(:),   INTENT(inout) :: xout
      REAL(wp), DIMENSION(:),   INTENT(inout) :: yout
      INTEGER,  DIMENSION(:),   INTENT(inout) :: idx
      INTEGER,                  INTENT(out)   :: nout

      INTEGER  :: nmax,  i, j, nseg, hnx, icount, nx0
      INTEGER  :: i0, i1, i2
      REAL(wp) :: yavg, ystd, val1, val2

      nmax  = SIZE(xin)

      ! make sure nx is an odd number
      IF (mod(nx, 2) /= 1) THEN
        nx = nx - 1
      ENDIF

      nseg = nint(1.0_wp * nmax / nx)

      hnx = (nx - 1) / 2

      icount = 1
      DO i = 1, nseg
        i0 = i * nx - hnx
        i1 = i0 - hnx
        i2 = i0 + hnx
        IF (i2 > nmax) THEN
           i2 = nmax
        ENDIF

        nx0 = (i2 - i1) + 1

        IF (i1 < i2) THEN
           yavg = sum(yin(i1:i2) ) / nx0
           ystd = sqrt(sum((yin(i1:i2) - yavg)**2) / nx0)
        ELSE
           yavg = yin (i1)
           ystd = 0.0_wp
        ENDIF
        DO j = i1, i2
          val1 = yavg - threshold * yavg
          val2 = yavg + threshold * yavg
          IF (yin(j) >= val1 .and. yin(j) <= val2) THEN

             yout(icount) = yin(j)
             xout(icount) = xin(j)
             idx(icount) = j
             icount = icount + 1

          ENDIF
        ENDDO
      ENDDO
      nout = icount - 1

END SUBROUTINE remove_spikes