! 
! 2008.1.7  debuged   I. Tamagawa
!
program test
implicit none

  integer, parameter :: n=10
  real :: x1(n), x2(n)
  integer :: i, errnum
  
  do i=1, n
    read(*,*) x1(i)
  end do

  call interpol( x1, x2, n, errnum )

  write(*,*) 'number of error data  = ', errnum
  do i=1, n
  	write(*,*) x1(i), x2(i)
  end do
end






logical function ok( x, count )
implicit none
  real :: x
  integer :: count

  if( x < -9900 ) then
  	ok = (5<0)
	count = count + 1
  else
        ok = (5>0)
  endif
end function


subroutine interpol( x1, x2, n, errnum )
implicit none
 real, intent(in)   :: x1(n)
 real, intent(out)  :: x2(n)
 integer, intent(in) :: n
 integer, intent(out) :: errnum

!  x1  ==> x2   (error data are linearly interpoleted)

 integer :: i, err1st, errlast, i1, i2
 integer :: count
 logical, external :: ok

 count = 0
 errnum = 0

 do i=1, n
	if( ok(x1(i), count) )  then
	  x2(i) = x1(i)
	  if( count > 0 ) then
	       errlast = i - 1
	       call linear_interpol( x1, x2, err1st, errlast, n )
	       errnum = errnum + (errlast-err1st+1)
	       count = 0
	  end if
	else
	  if( count == 1 ) then
	  	err1st = i
	  end if
	end if
 end do
 if( count > 0 ) then
    errlast = n
    call linear_interpol( x1, x2, err1st, errlast, n )
    errnum = errnum + count
 end if

end subroutine



subroutine find_next_ok( x, n, errlast, i2 )
! find next good data in x(i)   i > errlast+1

  implicit none
  integer, intent(in) :: n
  real, intent(in) :: x(n)
  integer :: errlast, i2
  integer :: i, dum
  logical, external :: ok

  !write(*,*) 'find next ok data'
  dum=0
  i = errlast+2
  do while( (.NOT. ok(x(i), dum)) .AND. (i<=n) )
  !write(*,*) 'i in while loop', i 
   i=i+1
  end do

  if( i == n ) then
    if( ok(x(i), dum) ) then
    	i2 = i
    else
        i2 = errlast + 1
    end if

  else
      i2 = i
  end if

  !write(*,*) ' i2 = ', i2
end


subroutine find_last_ok( x, n, err1st, i1 )
! find last good data in x(i)   i < err1st-1


  implicit none
  integer, intent(in) :: n
  real, intent(in) :: x(n)
  integer :: err1st, i1
  integer :: i, dum
  logical, external :: ok

  !write(*,*) 'find before  ok data'
  dum=0
  i = err1st-2
  do while( (.NOT. ok(x(i), dum)) .AND. (i>0) )
  !write(*,*) 'i in while loop', i 
   i=i-1
  end do

  if( i <= 0  ) then
        i1 = err1st - 1
  else
        i1 = i
  end if

  !write(*,*) ' i1 = ', i1
end











subroutine linear_interpol( x1, x2, err1st, errlast, n )
implicit none
  integer, intent(in) :: n
  real, intent(in) :: x1(n)
  real             :: x2(n)
  integer          :: err1st, errlast
  integer          :: i1, i2, i

  real ::  a, b


! write(*,*) ' interpol (err1st and errlast)', err1st, errlast

  if( err1st == 1 .AND. errlast == n ) then    ! NO GOOD DATA
    write(0,*) 'Give up to interporate error data, because all', &
      'data are error data'
    do i=1,n
     x2(i)=x1(i)
    end do

    return
  end if

  if( err1st > 1 .AND. errlast < n) then
     i1 = err1st - 1;  i2 =errlast +1

  else if( err1st == 1 .AND. errlast < n-1 ) then
        i1 = errlast +1 
	call find_next_ok( x1, n, errlast, i2 )
	errlast = i2 - 1


  else if( errlast == n .AND. err1st > 2 ) then
	i2 = err1st -1  	
  	call find_last_ok( x1, n, err1st, i1 )
	err1st = i1+1
  end if

  if( i2 > i1 ) then
    a = (x1(i2)-x1(i1))/(i2-i1)
    do i=err1st, errlast
	   x2(i) =  x1(i1) + a*(i-i1)
    end do
  else		!  only one datum is ok (maybe)
    write(0,*) 'Only one datum is OK.  Number of error data is wrong!'
    do i=1, n
    	x2(i) = x1(i1)
    end do
  end if
    	   

end subroutine



    
