c     use dflib
      character*5 dummy
      integer*4 h
      real*4, allocatable, dimension(:) :: x,y,m,q25,q75,wt,y_temp,
     &                                     wt_temp,y_star
      integer*4, allocatable, dimension(:) :: iperm
      open(2,file='clean.dat',status='old')
      open(3,file='clean.res',status='unknown')
      open(1,file='sub.res',status='unknown')
      n_data=0
      read(2,*) dummy
5     continue
         read(2,*,end=10) x_dummy
         n_data=n_data+1
         goto 5
10    continue
      write(6,*) 'ndata is:',n_data
      rewind(2)
      h=500
      allocate( y(n_data),x(n_data),m(n_data),q25(n_data),q75(n_data),
     &          wt(-h:h),y_temp(2*h+1),wt_temp(2*h+1),iperm(2*h+1),
     &          y_star(n_data) )
      read(2,*) dummy
      do i=1,n_data  
         read(2,*) x(i),y(i)
         end do
c
c     set up weights
c
      wt_tot=0.0
      do k=-h,h
         wt(k)=1.0 - (real(k)/real(h))**2
         wt_tot=wt_tot+wt(k)
         end do
      wt=wt/wt_tot
c
c     calculate quantiles
c
      do i=1,n_data
         i_low=max(1,i-h)
         i_up=min(n_data,i+h)
         j_max=i_up-i_low+1
         y_temp(1:j_max)=y(i_low:i_up)
         wt_temp(1:j_max)=wt(i_low-i:i_up-i)
         if ((i_low==1).or.(i_up==n_data)) then
            wt_tot=0.0
            do j=1,j_max
               wt_tot=wt_tot+wt_temp(j)
               iperm(j)=j
               end do
            wt_temp=wt_temp/wt_tot
         else
            do j=1,j_max
               iperm(j)=j
               end do
         end if
         call svrgp(j_max,y_temp(1:j_max),y_temp(1:j_max),
     &              iperm(1:j_max))
         cum_wt=0.0
         j=1
         do while (cum_wt<0.25)
            delta_w=wt_temp(iperm(j))
            cum_wt=cum_wt+delta_w
            j=j+1
            end do
         phi=(cum_wt-0.25)/delta_w
         q25(i)=phi*y_temp(j-2)+(1.0-phi)*y_temp(j-1)
         do while (cum_wt<0.5)
            delta_w=wt_temp(iperm(j))
            cum_wt=cum_wt+delta_w
            j=j+1
            end do
         phi=(cum_wt-0.5)/delta_w
         m(i)=phi*y_temp(j-2)+(1.0-phi)*y_temp(j-1)
         do while (cum_wt<0.75)
            delta_w=wt_temp(iperm(j))
            cum_wt=cum_wt+delta_w
            j=j+1
            end do
         phi=(cum_wt-0.75)/delta_w
         q75(i)=phi*y_temp(j-2)+(1.0-phi)*y_temp(j-1)
         y_star(i)=(y(i)-m(i))/(q75(i)-q25(i))
         write(3,*) x(i),y_star(i)
         write(1,*) x(i),y(i)-m(i)
         end do
100   format(1x,20000(f12.6,2x) )      
c
c     calculate interquartile range of standardized data
c
      call svrgn(n_data,y_star,y_star)
c
c     calculate variance of 1st normal using negative values
c
      xy=0.0
      ybar=0.0
      y2bar=0.0
      do i=1,n_data
         if (y_star(i)<0.0) then
            xy=xy+1.0
            ybar=ybar+(y_star(i)-ybar)/xy
            y2bar=y2bar+(y_star(i)**2-y2bar)/xy
         end if
         end do
      sigma=y2bar - ybar*ybar
      sigma=sqrt( sigma )/0.583819
      write(6,*) 'standard deviation of noise scale is:',sigma
      write(6,*) 's.e. via y2bar is:',sqrt(y2bar)
      sigma=sqrt(y2bar)
c      n_nonzero=j         
c      n_nonzero=n_data
c      i_q45=int( 0.45*real(n_nonzero) )
c      i_q55=int( 0.55*real(n_nonzero) )
c      q=(y_star(i_q55)-y_star(i_q45))/(anorin(0.55)-anorin(0.45))
c      write(6,*) 'q from 55-45 is:',q
c      i_q25=int( 0.25*real(n_nonzero) )
c      i_q75=int( 0.75*real(n_nonzero) )
c      q=(y_star(i_q75)-y_star(i_q25))/(anorin(0.75)-anorin(0.25))
c      write(6,*) 'q from 75 - 25 is:',q
      i_3q=0
      i_4q=0
      i_5q=0
      i_6q=0
      i_7q=0
      close(2)
      close(3)
      open(1,file='clean3.res',status='unknown')
      open(2,file='clean4.res',status='unknown')
      open(3,file='clean5.res',status='unknown')
      open(4,file='clean6.res',status='unknown')
      open(7,file='clean7.res',status='unknown')
      open(8,file='qqplot.res',status='unknown')
      do i=1,n_data
         if (y_star(i)>3.0*sigma) i_3q=i_3q+1
         if (y_star(i)>4.0*sigma) i_4q=i_4q+1
         if (y_star(i)>5.0*sigma) i_5q=i_5q+1
         if (y_star(i)>6.0*sigma) i_6q=i_6q+1
         if (y_star(i)>7.0*sigma) i_7q=i_7q+1
         ystar=(y(i)-m(i))/(q75(i)-q25(i))
         if (abs(ystar)>3.0*sigma) then
            write(1,*) x(i), ystar
         else
            write(1,*) x(i), 0.0d0
         end if
         if (abs(ystar)>4.0*sigma) then
            write(2,*) x(i), ystar
         else
            write(2,*) x(i), 0.0d0
         end if
         if (abs(ystar)>5.0*sigma) then
            write(3,*) x(i), ystar
         else
            write(3,*) x(i), 0.0d0
         end if
         if (abs(ystar)>6.0*sigma) then
            write(4,*) x(i), ystar
         else
            write(4,*) x(i), 0.0d0
         end if
         if (abs(ystar)>7.0*sigma) then
            write(7,*) x(i), ystar
         else
            write(7,*) x(i), 0.0d0
         end if
         write(8,*) anorin( real(i)/real(n_data+1) ),y_star(i)
         end do
      write(6,*) 'exceedances: 3sigma,4sigma,5sigma,6sigma,7sigma'
      write(6,*) i_3q,i_4q,i_5q,i_6q,i_7q
c     call setexitqq(QWIN$EXITNOPERSIST)
      end