subroutine addh(lheap,rface,nheap,ifnew)
      implicit double precision (a-h,o-z)
      dimension lheap(*),rface(*)
*
*     add to heap
*
      nheap=nheap+1
      lheap(nheap)=ifnew
      ipson=nheap
10    continue
      ipfath=ipson/2
      ifson=lheap(ipson)
      iffath=lheap(ipfath)
      if(rface(ifson).ge.rface(iffath)) goto 9000
      lheap(ipfath)=ifson
      lheap(ipson)=iffath
      ipson=ipfath
      if(ipson.gt.1) goto 10
9000  continue
      return
      end
      subroutine remh(lheap,rface,nheap,ifout)
      implicit double precision (a-h,o-z)
      dimension lheap(*),rface(*)
*
*     remove from heap
*
      ifout=lheap(1)
      lheap(1)=lheap(nheap)
      nheap=nheap-1
      ipfath=1
      ipson1=ipfath*2
10    continue
      ipson2=ipson1+1
      ifson1=lheap(ipson1)
      ifson2=lheap(ipson2)
      iffath=lheap(ipfath)
*
      rfath=rface(iffath)
      rson1=rface(ifson1)
      rson2=rface(ifson2)
*
      if(rfath.le.rson1.and.
     x rfath.le.rson2) ipexch=0
*
      if(rfath.ge.rson1.and.
     x rson1.ge.rson2) ipexch=ipson2
*
      if(rfath.ge.rson2.and.
     x rson2.ge.rson1) ipexch=ipson1
*
      if(rson1.ge.rfath.and.
     x rfath.ge.rson2) ipexch=ipson2
*
      if(rson2.ge.rfath.and.
     x rfath.ge.rson1) ipexch=ipson1
*
      if(ipexch.eq.0) goto 9000
      ifexch=ifson1
      if(ipexch.eq.ipson2) ifexch=ifson2
      lheap(ipexch)=iffath
      lheap(ipfath)=ifexch
      ipfath=ipexch
      ipson1=2*ipfath
      if(ipson1.le.nheap) goto 10
9000  continue
      return
      end