subroutine gvsort(r,ip,alow,ahi,istart,nentry,jp,n,nlim,toll)             
c...                                                                            
c... sort r (indexed by ip) by increasing value - vrs may 1988                  
c...                                                                            
      implicit double precision (a-h,o-z)
      dimension r(1),ip(1),alow(1),ahi(1),istart(1),nentry(1),jp(1)             
      al=r(1)                                                                   
      ip(1)=1                                                                   
      ah=al                                                                     
      do 1 l=2,n                                                                
      ip(l)=l                                                                   
      rrr=r(l)                                                                  
      al=min(al,rrr)                                                            
1     ah=max(ah,rrr)                                                            
      nseg=1                                                                    
      istart(1)=0                                                               
      nentry(1)=n                                                               
      alow(1)=al                                                                
      ahi(1)=ah                                                                 
999   al=alow(nseg)                                                             
      ah=ahi(nseg)                                                              
      if((ah-al).lt.toll)goto 888                                               
      nen=nentry(nseg)                                                          
      ist=istart(nseg)                                                          
      if(nen.lt.9)goto 777                                                      
c... bisect the list                                                            
      sss=(al+ah)*0.5d0                                                         
      ahi(nseg+1)=ah                                                            
      al2=ah                                                                    
      ah=al                                                                     
      ne=0                                                                      
      nep=0                                                                     
      do 5 l=1,nen                                                              
      k=ip(ist+l)                                                               
      rrr=r(k)                                                                  
      if(rrr.lt.sss)goto 6                                                      
      nep=nep+1                                                                 
      jp(nep)=k                                                                 
      al2=min(al2,rrr)                                                          
      goto 5                                                                    
6     ne=ne+1                                                                   
      ip(ist+ne)=k                                                              
      ah=max(ah,rrr)                                                            
5     continue                                                                  
      ahi(nseg)=ah                                                              
      nentry(nseg)=ne                                                           
      ist=ist+ne                                                                
      if(ist.ge.nlim)goto 999                                                   
      nseg=nseg+1                                                               
      istart(nseg)=ist                                                          
      nentry(nseg)=nep                                                          
      alow(nseg)=al2                                                            
      do 4 l=1,nep                                                              
4     ip(l+ist)=jp(l)                                                           
      goto 999                                                                  
c... short vector sort                                                          
777   ne=nen-1                                                                  
      do 2 l=1,ne                                                               
      k=ip(l+ist)                                                               
      rrr=r(k)                                                                  
      nep=l+1                                                                   
      do 3 m=nep,nen                                                            
      j=ip(m+ist)                                                               
      sss=r(j)                                                                  
      if(sss.ge.rrr)goto 3                                                      
      ip(m+ist)=k                                                               
      k=j                                                                       
      rrr=sss                                                                   
3     continue                                                                  
2     ip(l+ist)=k                                                               
888   nseg=nseg-1                                                               
      if(nseg.ne.0)goto 999                                                     
      return                                                                    
      end