FUNCTION NCTREE(X,DX, INX,NX,VECT,NDIM,NVEC,INDLR)
*     binary search tree - returns number NC, and the indices in common
*
*     A search is made for key vectors, which satisfy the conditions.
*     The arguments INX,NX, VECT,NDIM,NVEC, INDLR are as before in DFTREE.
*     The user has to specify the limits in the arrays X(NDIM),DX(NDIM).
*
*     The returned function value NC is the number of key vectors found,
*     satisfying the conditions given by the limits.
*     The range search results of the function NCTREE is stored in the
*     integer array INDC(.) in the common /CINDEX/. The default dimension
*     of the array INDC(.) is 100, sufficient for up to 100 indices found
*     in a range search. The user may define a larger INDC(.) array dimension
*     NCDIM in the common, which has to be specified as the last argument
*     in the call.
*     The indices of the found key vectors are in the array
*                        INDC(J), J = 1, NC
*     in the common /CINDC/.
*     Additional information on the range search in the common is:
*        NTRY   = number of test attempts (measures the efficiency)
*        ISTMAX = maximum length of the stack, used in the range search
*     The default dimension of the stack is 100, which should be sufficient
*     for almost all cases. The programs stops in case of a stack overflow.
*
      REAL X(*),DX(*),VECT(NDIM,NVEC)
      INTEGER INX(NX),INDLR(2,NVEC)         ! index vector, tree pointer
      PARAMETER               (MINDC=100)   ! default dimension of
      COMMON/CINDEX/NINDC,NTRY,ISTMAX,INDC(MINDC)  ! result index common
      PARAMETER  (NSTACK=100)
      INTEGER KST(NSTACK),LST(NSTACK)       ! stack
      LOGICAL TXL,TXR
*     ...
      NTRY=0
      NC=0
      IST=1
      KST(IST)=1                     ! push
      LST(IST)=1                     ! 1. component
      ISTMAX=0
  10  I=KST(IST)                     ! pop element
      L=LST(IST)                     ! index
      IST=IST-1
      IF(I.EQ.0) STOP 'I is zero '
      TXL=   VECT(INX(L),I).GE.X(L)-DX(L)
      TXR=   VECT(INX(L),I).LE.X(L)+DX(L)
      NTRY=NTRY+1                    ! count number of attempts
      IF(TXL.AND.TXR) THEN
         DO K=1,NX                     ! check rectangle
          IF(ABS(VECT(INX(K),I)-X(K)).GT.DX(K)) GOTO 20
         END DO
         NC=NC+1                        ! count
         IF(NC.LE.NINDC) INDC(NC)=I     ! save index
      END IF
  20  IF(INDLR(2,I).NE.0.AND.TXR) THEN
         IST=IST+1
         IF(IST.GT.NSTACK) STOP 'NCTREE stack overflow'
         KST(IST)=INDLR(2,I)         ! push r
         LST(IST)=MOD(L,NX)+1        ! next component
         ISTMAX=MAX(ISTMAX,IST)
      END IF
      IF(INDLR(1,I).NE.0.AND.TXL) THEN
         IST=IST+1
         IF(IST.GT.NSTACK) STOP 'NCTREE stack overflow'
         KST(IST)=INDLR(1,I)         ! push l
         LST(IST)=MOD(L,NX)+1        ! next component
         ISTMAX=MAX(ISTMAX,IST)
      END IF
      IF(IST.NE.0) GOTO 10
      NCTREE=NC
      END