C     ******************************************************************
C     ******************************************************************

      subroutine iniset()

      implicit none

C     This subroutine initiates the set S of subboxes.

C     PARAMETERS
      include "dim.f"

C     COMMON SCALARS
      integer qf

C     COMMON ARRAYS
      double precision subboxlb(qmax),subboxl(lmax),subboxu(lmax)

C     COMMON BLOCKS
      common /setS/ subboxlb,subboxl,subboxu,qf

      qf = 0

      end 

C     ******************************************************************
C     ******************************************************************

      subroutine putsubbox(n,l,u,lb,overflow)

      implicit none

C     SCALAR ARGUMENTS
      integer n
      logical overflow
      double precision lb

C     ARRAY ARGUMENTS
      double precision l(n),u(n)

C     This subroutine add a subbox to S.

C     PARAMETERS
      include "dim.f"

C     COMMON SCALARS
      integer qf

C     COMMON ARRAYS
      double precision subboxlb(qmax),subboxl(lmax),subboxu(lmax)

C     LOCAL SCALARS
      integer i,ind,ind1,ind2,j,k

C     COMMON BLOCKS
      common /setS/ subboxlb,subboxl,subboxu,qf

      if ( qf + 1 .gt. qmax .or. ( qf + 1 ) * n .gt. lmax ) then
          overflow = .true.
          return
      end if

C     Add the new subbox at the end

      qf = qf + 1

      subboxlb(qf)  = lb

      do i = 1,n
          ind = ( qf - 1 ) * n + i
          subboxl(ind) = l(i)
          subboxu(ind) = u(i)
      end do

      overflow = .false.

C     Move the new last subbox up

      j = qf

 10   continue

      k = j / 2

      if ( k .ge. 1 ) then
          if ( subboxlb(j) .lt. subboxlb(k) ) then
              call ddswap(subboxlb(j),subboxlb(k))

              do i = 1,n
                  ind1 = ( j - 1 ) * n + i
                  ind2 = ( k - 1 ) * n + i
                  call ddswap(subboxl(ind1),subboxl(ind2))
                  call ddswap(subboxu(ind1),subboxu(ind2))
              end do

              j = k
              go to 10
          end if
      end if

      end

C     ******************************************************************
C     ******************************************************************

      subroutine getsubbox(n,l,u,lb,underflow)

      implicit none

C     SCALAR ARGUMENTS
      integer n
      logical underflow
      double precision lb

C     ARRAY ARGUMENTS
      double precision l(n),u(n)

C     This subroutine remove a subbox from S.

C     PARAMETERS
      include "dim.f"
 
C     COMMON SCALARS
      integer qf

C     COMMON ARRAYS
      double precision subboxlb(qmax),subboxl(lmax),subboxu(lmax)

C     LOCAL SCALARS
      integer i,ind,ind1,ind2,j,k

C     COMMON BLOCKS
      common /setS/ subboxlb,subboxl,subboxu,qf

      if ( qf .le. 0 ) then
          underflow = .true.
          return
      end if

C     Remove the first subbox and put the last one in its place

      lb  = subboxlb(1)

      do i = 1,n
          l(i) = subboxl(i)
          u(i) = subboxu(i)
      end do

      do i = 1,n
          ind = ( qf - 1 ) * n + i
          subboxl(i) = subboxl(ind)
          subboxu(i) = subboxu(ind)
      end do

      subboxlb(1)  = subboxlb(qf)

      qf = qf - 1

      underflow = .false.

C     Move down the new first subbox

      j = 1

 10   continue

      k = 2 * j

      if ( k .le. qf ) then
          if ( k .lt. qf ) then
              if ( subboxlb(k + 1) .lt. subboxlb(k) ) then
                  k = k + 1
              end if
          end if

          if ( subboxlb(j) .gt. subboxlb(k) ) then
              call ddswap(subboxlb(j),subboxlb(k))

              do i = 1,n
                  ind1 = ( j - 1 ) * n + i
                  ind2 = ( k - 1 ) * n + i
                  call ddswap(subboxl(ind1),subboxl(ind2))
                  call ddswap(subboxu(ind1),subboxu(ind2))
              end do

              j = k
              go to 10
          end if
      end if

      end

C     ******************************************************************
C     ******************************************************************

      subroutine ddswap(a,b)

      implicit none

C     SCALAR ARGUMENTS
      double precision a,b

C     LOCAL SCALARS
      double precision c

      c = a
      a = b
      b = c

      end