!Identifying the flat portion of a curve by computing the mean slope of
!all possible pairs. Data must be ranked with xvar from low to high
      subroutine findplateau(npoints,xvar,yvar,iflat,
     &  flatmean,radius)
      implicit none
      integer npoints,iflat(npoints),isoutlier_1side
      double precision xvar(npoints),yvar(npoints),
     &  flatmean,radius,t0,ymax
      integer i,j,k,m,n,nstart,nend,isslopezero,no1,no2,
     &  no3,nset
      double precision der(npoints*(npoints-1)/2),std,term,
     &  student_t,slopemin,Sign_Level,coeffvar,coeffvarmin,
     &  xx(npoints),yy(npoints),ainter,bslope,r,eps
      parameter(Sign_Level=0.90d0,eps=1.0d-6,nset=3)
      
      do i=1,npoints
        iflat(i)=1
      enddo
      nend=-9999
      nstart=0
      slopemin=1.0d+9
      coeffvarmin=1.0d+9
      flatmean=-9999.0d0
      radius=-9999.0d0
      
      ymax=0.0d0
      do i=1,nset
        ymax=ymax+yvar(i)
      enddo
      ymax=ymax/dble(nset)
      nstart=1
      do i=2,npoints-nset+1
        term=0.0d0
        do j=i,i+nset-1
          term=term+yvar(j)
        enddo
        term=term/dble(nset)
        if(term.ge.ymax)then
          ymax=term
          nstart=i
        endif
      enddo
      no1=nstart
      do i=nstart+1,nstart+nset-1
        if(yvar(i).ge.yvar(no1))then
          no1=i
        endif
      enddo
      nstart=1
500   k=no1-nstart+1
      if(k.lt.3)then
        nstart=no1
        goto 510
      endif
      do i=nstart,no1
        xx(i-nstart+1)=xvar(i)
        yy(i-nstart+1)=yvar(i)
      enddo
      call linregres(k,xx,yy,Sign_Level,ainter,
     &  bslope,r,radius,isslopezero)
      if(isslopezero.eq.-1)then
        nstart=nstart+1
        goto 500
      endif
510   nend=no1
      if((nend-nstart+1).le.4)then
        k=nstart
        do i=1,k-1
          do j=k,nend
            if(yvar(j).lt.yvar(i))then
              nstart=no1
            endif
          enddo
        enddo
      endif
      nend=npoints
600   k=nend-no1+1
      if(k.lt.3)then
        nend=no1
        goto 610
      endif
      do i=no1,nend
        xx(i-no1+1)=xvar(i)
        yy(i-no1+1)=yvar(i)
      enddo
      call linregres(k,xx,yy,Sign_Level,ainter,
     &  bslope,r,radius,isslopezero)
      if(isslopezero.eq.-1)then
        nend=nend-1
        goto 600
      endif
610   if((nend-no1+1).le.4)then
        k=nend
        do i=k+1,npoints
          do j=no1,k
            if(yvar(j).lt.yvar(i))then
              nend=no1
            endif
          enddo
        enddo
      endif
      if((nend-nstart).eq.0)then
        if(no1.eq.npoints)return
        if(no1.eq.1)then
          do i=1,npoints
            iflat(i)=3
          enddo
          return
        endif
        nstart=no1-1
        nend=no1
      endif  
      goto 100        

!find the three largest value in yvar
      if(yvar(1).lt.yvar(2))then
        no1=2
        no2=1
      else
        no1=1
        no2=2
      endif
      if(yvar(3).ge.yvar(no1))then
        no3=no2
        no2=no1
        no1=3
      else
        if(yvar(3).lt.yvar(no2))then
          no3=3
        else
          no3=no2
          no2=3
        endif
      endif
      do i=4,npoints
        if(yvar(i).ge.yvar(no1))then
          no3=no2
          no2=no1
          no1=i
        else
          if(yvar(i).ge.yvar(no2))then
            no3=no2
            no2=i
          else
            if(yvar(i).ge.yvar(no3))then
              no3=i
            endif
          endif
        endif
      enddo
      If(no1.gt.no2.and.no1.gt.no3)then
        nend=no1
        if(no2.gt.no3)then
          nstart=no3
        else
          nstart=no2
        endif
      endif
      If(no2.gt.no1.and.no2.gt.no3)then
        nend=no2
        if(no1.gt.no3)then
          nstart=no3
        else
          nstart=no1
        endif
      endif
      If(no3.gt.no1.and.no3.gt.no2)then
        nend=no3
        if(no1.gt.no2)then
          nstart=no2
        else
          nstart=no1
        endif
      endif
      if(nstart.gt.1)then
        i=1
2       if(dabs(yvar(i)-yvar(no1)).lt.eps.or.
     &    dabs(yvar(i)-yvar(no2)).lt.eps.or.
     &    dabs(yvar(i)-yvar(no3)).lt.eps)then
          nstart=i
          goto 1
        endif
        i=i+1
        if(i.eq.nstart)goto 1
        goto 2
      endif
1     if(nend.lt.npoints)then
        i=nend+1
3       if(dabs(yvar(i)-yvar(no1)).lt.eps.or.
     &    dabs(yvar(i)-yvar(no2)).lt.eps.or.
     &    dabs(yvar(i)-yvar(no3)).lt.eps)then
          nend=i
        endif
        if(i.eq.npoints)goto 4
        i=i+1
        goto 3
      endif
4     do i=nstart,nend
        xx(i-nstart+1)=xvar(i)
        yy(i-nstart+1)=yvar(i)
      enddo
      k=nend-nstart+1
      call linregres(k,xx,yy,Sign_Level,ainter,
     &  bslope,r,radius,isslopezero)
      if(isslopezero.eq.-1)then
        if(no1.eq.npoints)return
        if(no1.eq.nstart)then
          nstart=no1
          nend=no1
          goto 100
        endif
        if(no1.gt.no2)then
          nstart=no2
          nend=no1
        else
          nstart=no1
          nend=no2
        endif
        goto 100
      endif
      slopemin=dabs(bslope)
      call stdmean(k,yy,std,flatmean)
      coeffvarmin=std/flatmean
10    if(nend.eq.npoints)goto 50
      if(no1.eq.nend)then
        if(dabs(yvar(nend+1)-yvar(no1)).lt.eps.or.
     &    dabs(yvar(nend+1)-yvar(no2)).lt.eps.or.
     &    dabs(yvar(nend+1)-yvar(no3)).lt.eps)then
          nend=nend+1
          goto 10
        else
          goto 50
        endif
      endif
      do k=no1,nend+1
        xx(k-no1+1)=xvar(k)
        yy(k-no1+1)=yvar(k)
      enddo
      k=nend+1-no1+1
      call linregres(k,xx,yy,Sign_Level,ainter,bslope,
     &  r,radius,isslopezero)
      if(isslopezero.eq.-1)goto 50
      call stdmean(k,yy,std,flatmean)
      coeffvar=std/flatmean
      nend=nend+1
      goto 10
50    if(nstart.eq.1)goto 100
      if(no1.eq.nstart)then
        if(dabs(yvar(nstart-1)-yvar(no1)).lt.eps.or.
     &    dabs(yvar(nstart-1)-yvar(no2)).lt.eps.or.
     &    dabs(yvar(nstart-1)-yvar(no3)).lt.eps)then
          nstart=nstart-1
          goto 50
        else
          goto 100
        endif
      endif
      do k=nstart-1,no1
        xx(k-nstart+1+1)=xvar(k)
        yy(k-nstart+1+1)=yvar(k)
      enddo
      k=no1-(nstart-1)+1
      call linregres(k,xx,yy,Sign_Level,ainter,bslope,
     &  r,radius,isslopezero)
      if(isslopezero.eq.-1)goto 100
      call stdmean(k,yy,std,flatmean)
      coeffvar=std/flatmean 
      nstart=nstart-1
      goto 50

100   if(nstart.le.nend)then
        n=0
        do i=nstart,nend
          n=n+1
          der(n)=yvar(i)
        enddo
        if(n.gt.1)then
          call stdmean(n,der,std,flatmean)
          t0=student_t(n-1,Sign_Level)
          radius=t0*std/dsqrt(dble(n))
        else
          flatmean=der(1)
          radius=-9999.0d0
        endif
        do i=1,nstart-1
          iflat(i)=1
        enddo
        do i=nstart,nend
          iflat(i)=2
        enddo
        do i=nend+1,npoints
          iflat(i)=3
        enddo
        if(n.gt.2)then
          k=isoutlier_1side(n,der,-1)
          if(k.eq.1)then
            iflat(nstart)=1
          endif
          if(k.eq.n)then
            iflat(nend)=3
          endif
        endif
      endif
      return
      end