subroutine gvlst2(alat,plat,q,n1,n2,n3,gmin,gmax,mshlst,job,ngmx, ng,kv,gv,igv)!Setup a list of recip vectors within cutoff |q+G| < gmax
  use m_ftox
  use m_shortn3,only: shortn3_initialize,shortn3,nout,nlatout
  use m_lgunit,only:stdo
  use m_mksym,only:  ngrp=>lat_nsgrp,gsym=>rv_a_osymgr
  implicit none
  intent(in)::    alat,plat,q,n1,n2,n3,gmin,gmax,mshlst,job,ngmx
  intent(out)::                                                      kv,gv,igv 
  intent(inout)::                                                 ng
  !i   alat     Lattice constant
  !i   plat     Real-space primitive lattice vectors
  !i   q        vectors |q+G|<gmax are included in list
  !i   mshlst   if first entry is nonzero, a list of allowed values
  !i            n1,n2,n3 may take.  First entry is the size of
  !i            the list; then follows the list itself.
  !i   job      1s digit
  !i              0 return ng only
  !i              1 return kv and igv
  !i              8 return kv and gv, and then sorted
  !i                any combination of 1,8 is allowed
  !i          +500: gv without q 
  !i          +1000: symcheck
  !i   ngmx     Leading dimension of kv,gv,igv
  !i   gmin     Lower bound for reciprocal lattice vectors, in a.u.
  !i  n1..3    On input, max # divisions along the three lattice vectors.
  !i           (It is an error for both gmax and n1..n3 to be zero.)
  !i           Otherwise, input n1..n3 additionally constrain which
  !i           vectors are added to list; see Remarks.
  !i
  !i   gmax    On input, cutoff for reciprocal lattice vectors, in a.u.
  !i           Energy cutoff is gmax**2.
  !i           If input gmax is zero, gvlst2 will generate it from n1..n3
  !i           (It is an error for both gmax and n1..n3 to be zero.)
  !io   ng       Number of lattice vectors
  !o   kv       indices for gather/scatter operations.
  !o            kv(ig,i=1,2,3) for vector ig point to which entry
  !o            (i1,i2,i3) vector ig belongs
  !o   gv       list of reciprocal lattice vectors G
  !o   igv      list of reciprocal lattice vectors G, represented as
  !o            three integers (the multiples of qlat)
  !o            gv and igv are related by:
  !o              gv(1:3,1:ng) = 2*pi/alat * (qlat * igv(1:ng))
  !r Remarks
  !r   Collects a list of q + reciprocal lattice vectors (G+q) that lie
  !r   within a cutoff gmax.  List is optionally sorted by length (logical: lsort)
  !r   Vectors G are integer multiples of primitive r.l.v. qlat.
  !r
  !r   Additional constraints may be imposed, namely that the number of
  !r   multiples in each of axes 1..3 not exceed input values n1..n3.
  !r
  !r   For copying the list of points to/from a regular mesh, use:
  !r     call gvgetf(ng,1,kv,k1,k2,k3,c,c0)
  !r       to collect elements into list c0 from 3D array c
  !r     call gvputf(ng,1,kv,k1,k2,k3,c0,c)
  !r       to poke elements from list c0 into 3D array c
  !r     call gvaddf(ng,kv,k1,k2,k3,c0,c)
  !r       to add elements from list c0 into 3D array c
  !r
  !u Updates See git log after 2009
  !u   11 Jul 08 New argument gmin
  !u   01 Jun 01 revised gvlst2 together with gvlist.  They form similar
  !u             operations but with different functions; see gvlist.f
  !u   26 Mar 01 Another bug fix for input n1..n3 and gmax nonzero
  !u   06 Mar 01 Bug fix for input n1..n3 and gmax nonzero
  !u   Routine was adapted from T. Kotani, routine getgv2
  ! ----------------------------------------------------------------------
  integer :: n1,n2,n3,ng,ngmx,kv(ngmx,3),igv(ngmx,3),mshlst(0:*) !,igv2(3,*)
  double precision :: alat,gmin,gmax,gv(ngmx,3),plat(3,3),q(3)
  integer :: ig,n1max,n1min,n2max,n2min,n3max,n3min,i1,i2,i3,nn,i
  integer :: n1l,n2l,n3l
  integer :: PRTG,PRTG2,iset(3),ipr,job,job0,job1,job8,k1,k2,k3
  double precision :: qlat(3,3),vol,pi,tpiba,qpg(3),q2
  double precision :: gmin0,gmax0,gmin2,gmax2,h1,h2,h3,ddot,tol
  character(256) :: outs
  parameter (PRTG=30,PRTG2=100,tol=1d-8)
  real(8),parameter:: tolg2=1d-2
  logical::  lgpq,lgv,lsort,ligv !,ligv2
  real(8):: plat1(3,3),qlat1(3,3),gg,gs(3)
  integer:: j1,j2,j3,m,jj1,jj2,jj3,nn1,nn2,nn3,i123(3),jjj(3),jg,igrp,jjg
  real(8):: rlatp(3,3),xmx2(3),gvv(3),diffmin
  integer :: nginit,kv_tmp(ngmx,3),igv_tmp(ngmx,3),ips(ngmx),jx,nxx,itemp(ngmx),ix,iprint
  real(8):: gv_tmp(ngmx,3)
  call tcn('gvlst2')
  call pshpr(iprint()-30)
  call getpr(ipr)
  call dinv33(plat,1,qlat,vol)
  pi = 4d0*datan(1d0)
  tpiba = 2*pi/alat
  gmin0  = gmin/tpiba
  gmax0  = gmax/tpiba
  if (gmin < 0)  call rx('gvlst2: input gmin <= 0')
  if (gmax <= 0) call rx('gvlst2: input gmax <= 0')
  job0 = mod(job,100)
  ligv  = mod(job0,2)/=0    !1,9 !Get igv
  lgv   = mod(job0/4,4)/=0  !8,9 !Get gv
  lsort = mod(job0/8,2)/=0  !8,9 sorted
  lgpq  = mod(job/100,10)>4 !+500 or not
  
  !! ... Basis vectors for real-space mesh and recip-space supercell
  nn1=n1
  nn2=n2
  nn3=n3
  n1min=0; n1max=nn1-1
  n2min=0; n2max=nn2-1
  n3min=0; n3max=nn3-1
  if(nn1*nn2*nn3==0) then
     print *,'nn1,nn2,nn3=',nn1,nn2,nn3
     call rx('gvlst2: require n1*n2*n3/=0')
  endif   
!  if(nn1==0) call gvlstn(qlat(1,1),qlat(1,2),qlat(1,3),q,mshlst,gmax0,nn1) 
!  if(nn2==0) call gvlstn(qlat(1,2),qlat(1,3),qlat(1,1),q,mshlst,gmax0,nn2) 
!  if(nn3==0) call gvlstn(qlat(1,3),qlat(1,1),qlat(1,2),q,mshlst,gmax0,nn3) 
  do m = 1, 3
     plat1(m,:) = plat(m,:)*[1d0/nn1,1d0/nn2,1d0/nn3]
     qlat1(m,:) = qlat(m,:)*[nn1,nn2,nn3]
  enddo
  n1l=nn1
  n2l=nn2
  n3l=nn3
  gmax2 = (gmax0-tol)**2 !this tol is needed for gmax given by gvctof
  gmin2 = gmin0**2
  call shortn3_initialize(qlat1) !initialization for m_shoten3 for qlat1
  ig=0
  j1loop: do j1 = 0,nn1-1 !! --- g vectors, shorten, count and keep if within gmax ---
     j2loop: do j2 = 0,nn2-1 
        j3loop: do j3 = 0,nn3-1
           jjj=[j1,j2,j3] 
           qpg= [1d0/nn1, 1d0/nn2, 1d0/nn3] * (jjj+ matmul(q,plat(:,:))) !frac corrdinate on qlat1.
           call shortn3(qpg) ! return nout,nlatout
           gs = matmul(qlat1(:,:), (qpg+nlatout(:,1)))
           gg = sum(gs**2)
           if(gmin2<= gg .AND. gg < gmax2) then
              ig = ig+1
              if (job0==0) cycle
              if (ig > ngmx) write(stdo,*) 'ERROR: ig ngmx=',ig,ngmx
              if (ig > ngmx) call rx('gvlist2: ng exceeds ngmx')
              kv_tmp(ig,:) = jjj+1
              if(ligv) igv_tmp(ig,:) = jjj + [nn1,nn2,nn3]*nlatout(:,1) !integer index of Gvec
              if(lgv)  gv_tmp(ig,:)= merge(gs,gs-q,mask=lgpq)
           endif
        enddo j3loop
     enddo j2loop
  enddo j1loop
  ng = ig
  if(job0/=0) then !symmetry checker
     nginit = ng
     ng=0
     ips=0
     do ig = 1,nginit
        if(job>999.and.ips(ig)==0) then 
           itemp=0
           ix=0
           do igrp = 1, ngrp
              gvv = matmul(gsym(:,:,igrp),gv_tmp(ig,:)) !  ... gvv = g(k) gv
              !write(6,ftox) 'iii ig igrp gvv=',ig,igrp,ftof(gvv)
              jg = findloc([(sum(abs(gvv-gv_tmp(jg,:)))<tolg2,jg=1,nginit)],value=.true.,dim=1)
              if(jg/=0) then
                 ix=ix+1
                 itemp(ix)=jg
              else
                 goto 70
              endif   
           enddo
           ips(itemp(1:ix))=1 !jg=itemp(1:ix) need to be included.
        endif
        ng = ng+1
        kv(ng,:) = kv_tmp(ig,:)
        if(ligv) igv(ng,:)=igv_tmp(ig,:)
        if(lgv )  gv(ng,:)= gv_tmp(ig,:) 
70      continue
     enddo !   write(stdo,ftox)'gmax ng nginit ngmx=',ftof(gmax),ng,nginit,ngmx,ftof(q)
  endif
  if(lsort) then
     gvsort:block !Sort vectors -- !call gvlsts(ng,gv(1:ng,1:3),kv(1:ng,1:3),igv(1:ng,1:3),ligv) 
       integer:: iprm(ng)
       call dvshel(1,ng, sum(gv(1:ng,1:3)**2,dim=2)*[((1d0 + 1d-15*ig),ig=1,ng)], iprm,1)
       gv(1:ng,1:3) = gv(iprm+1,1:3)
       kv(1:ng,1:3) = kv(iprm+1,1:3)
       if(ligv) igv(1:ng,1:3) =igv(iprm+1,1:3)
     endblock gvsort !write(stdo,ftox)'gmax gv=',ftof(gv(1,1:3))
  endif
  if(ipr>=PRTG)write(stdo,ftox)'gvlst2: gmax=',ftof(gmax,3),'a.u. created',ng,'vectors of',n1l*n2l*n3l,&
       '(',(ng*100)/(n1l*n2l*n3l),'%)'
  if(ipr >= PRTG2 .AND. ng > 0 .AND. ligv) then
     write(stdo,"(' G vectors (multiples of reciprocal lattice vectors)'/ '   ig    G1   G2   G3     E')")
     do  ig = 1, ng
        i123 = igv(ig,:)
        qpg= q + matmul(qlat,i123) 
        write(stdo,"(i5,1x,3i5,2x,f8.4)")  ig,i123,sum(qpg**2) *tpiba**2
        write(stdo,"('q  qpg=',3d13.5,3x,3d13.5)")  q,qpg
     enddo
  endif
  call poppr
  call tcx('gvlst2')
end subroutine gvlst2
subroutine gvgetf(ng,n,kv,k1,k2,k3,c,c0)!- Gathers Fourier coefficients from 3D array c into list c0.
  implicit none
  integer :: ng,n,k1,k2,k3,kv(ng,3)
  complex(8):: c0(ng,n),c(k1,k2,k3,n)
  integer :: ig,i,j1,j2,j3
  do i=1,n
     c0(:,i) = [(c(kv(ig,1),kv(ig,2),kv(ig,3),i), ig=1,ng)]
  enddo   
end subroutine gvgetf
subroutine gvputf(ng,n,kv,k1,k2,k3,c0,c)!- Pokes Fourier coefficients from list c0 into 3D array c.
  implicit none
  integer :: ng,n,k1,k2,k3,kv(ng,3)
  complex(8):: c0(ng,n),c(k1,k2,k3,n)
  integer :: ig,i,j1,j2,j3
  c=0d0
  do ig=1,ng
     c(kv(ig,1),kv(ig,2),kv(ig,3),:) = c0(ig,:)
  enddo   
end subroutine gvputf
subroutine gvaddf(ng,kv,k1,k2,k3,c0,c)! Adds Fourier coefficients from list c0 into 3D array c.
  implicit none
  integer :: ng,k1,k2,k3,kv(ng,3)
  complex(8):: c0(ng),c(k1,k2,k3)
  integer :: ig,j1,j2,j3
  do  ig = 1, ng
     j1 = kv(ig,1)
     j2 = kv(ig,2)
     j3 = kv(ig,3)
     c(j1,j2,j3) = c(j1,j2,j3) + c0(ig)
  enddo
end subroutine gvaddf
subroutine mshsiz(alat,plat,gmax,ngabc,ng) !Finds dimensions for a mesh of G vectors that satisfy a cutoff
  use m_lgunit,only:stdo
  !i Inputs
  !i   alat  :length scale of lattice and basis vectors, a.u.
  !i   plat  :primitive lattice vectors, in units of alat
  ! io  gmax  :On input, cutoff in reciprocal lattice vectors
  ! o        :Energy cutoff is gmax**2.
  ! o        :If input gmax is zero, it will be generated it from ngabc
  ! o        :(It is an error for both gmax and any of ngabc to be zero.)
  ! o        :On output, gmax is computed if input n1..n3 are nonzero,
  ! o        :or if job = 1
  ! io  ngabc :On input, max # divisions along the three recip lattice vecs.
  ! o        :(It is an error for both gmax and any of ngabc to be zero.)
  ! o        :Otherwise, those input ngabc initially zero are found that
  ! o        :satisfy cutoff criterion.
  ! o   ng    :number of G-vectors
  implicit none
  integer :: ngabc(3),ng,job,j
  double precision :: alat,plat(3,3),gmax
  logical :: change
  integer :: npfac,pfac(10),fmax,i,indx,iprint,i1,i2,i3,nn,nmin,nmx(3),nmxn(3),k,ngabcn(3),ngn,nginit,PRTG,kxx(1,1)
  double precision :: gmaxn,xx,q(3),gmax0,qlat(3,3),tpiba,facg,tolg,s1,s2,s3,ddot,gxx(1,1)
  character(256) :: outs
  parameter (fmax=600,tolg=1d0,PRTG=30)
  integer :: mshlst(0:3*fmax)
  logical:: fullmesh
  call tcn('mshsiz')
  q=0d0
  if (gmax == 0) then
     call gvctof(0,alat,plat,q,ngabc(1),ngabc(2),ngabc(3),gmax,ng)
     goto 999
  endif
  tpiba = 2*4d0*datan(1d0)/alat
  gmax0  = gmax/tpiba ! ... gmax0 = dimensionless gmax
  call dinv33(plat,1,qlat,xx)
  mshlst(0) = 0
  call ppfac(fmax,1,mshlst(1),mshlst)  !list of all allowed values of n1..n3
  if (mshlst(1) == 0) call rx('mshsiz: null set of allowed mesh points')
  ! --- Upper bound for n1..n3 : guaranteed to hold all G<gmax ---
  do  i = 1, 3
     i2 = mod(i,3)+1
     i3 = mod(i2,3)+1
     call gvlstn(qlat(1,i),qlat(1,i2),qlat(1,i3),q,mshlst,gmax0,nn)
     nmx(i)=nn
     if (ngabc(i) == 0) ngabc(i) = nn
  enddo
  gmaxn = gmax0
  facg = .995d0
  call pshpr(iprint()-30)
  call gvlst2(alat,plat,q,ngabc(1),ngabc(2),ngabc(3), 0d0,gmax,[0],000, 0,ng,kxx,gxx,kxx,kxx)
  ! ... Count the number of G vectors ng for initial n1..n3
  nginit = ng
  Getngloop: do !Reduce ng slowly until at least one of n1..n3 changes
     do  
        gmaxn = gmaxn * facg
        change = .false.
        do  i = 1, 3
           i2 = mod(i,3)+1
           i3 = mod(i2,3)+1
           call gvlstn(qlat(1,i),qlat(1,i2),qlat(1,i3),q,mshlst,gmaxn,nn)
           nmxn(i)=nn
           if (nmxn(i) /= nmx(i)) change = .TRUE. 
           !       The granularity of gvlstn may be too coarse.
           !       Don't assign, ngabcn(i) = nn but find next one smaller in mshlst
           indx = 1
           !call hunti(mshlst(1),mshlst,ngabc(i),0,indx)
           indx= findloc([(mshlst(j)<ngabc(i).and.ngabc(i)<=mshlst(j+1),j=1,mshlst(0)-1)],value=.true.,dim=1)
           indx = max(indx,1)
           ngabcn(i) = mshlst(indx)
        enddo
        if (change) exit
     enddo
     if(fullmesh()) then
        ngabc=ngabcn
        gmax=1d10
        ng=ngabc(1)*ngabc(2)*ngabc(3)
        exit
     endif   ! ... Count the number of G vectors for (smaller) trial ngabcn
     call gvlst2(alat,plat,q,ngabcn(1),ngabcn(2),ngabcn(3), 0d0,gmax,[0],000, 0,ngn,kxx,gxx,kxx,kxx)
     if (dble(ngn) >= nginit*tolg) then
        ng = ngn
        ngabc=ngabcn
     else
        exit
     endif
  enddo Getngloop 
  call poppr
  i1 = ngabc(1)
  i2 = ngabc(2)
  i3 = ngabc(3)
  if (max0(i1,i2,i3) == fmax .AND. iprint() >= 10) then
     write(stdo,301)
     write(stdo,300) fmax
     write(stdo,301)
  endif
300 format(' WARNING!'/' At least one of the mesh divisions reached its maximal value fmax = ',i4/ &
       ' You might need to increase parameter fmax in mshsiz')
301 format(/1x,79('*')/)
  if (iprint() >= PRTG) then
     s1 = alat*sqrt(ddot(3,plat(1,1),1,plat(1,1),1))/i1
     s2 = alat*sqrt(ddot(3,plat(1,2),1,plat(1,2),1))/i2
     s3 = alat*sqrt(ddot(3,plat(1,3),1,plat(1,3),1))/i3
     write(stdo,"('MSHSIZ: mesh has ',i0,' x ',i0,' x ',i0,' divisions; length =',3f10.3)")i1,i2,i3,s1,s2,s3
     write(stdo,"('      generated from gmax (a.u.)=',f12.4,': ',i0,' vectors of ', &
          i0,' (',i0,'%)')")  gmax,ng,i1*i2*i3,(ng*100)/(i1*i2*i3)
  endif
999 continue
  call tcx('mshsiz')
end subroutine mshsiz
subroutine gvlstn(q0,q1,q2,qp,mshlst,gmax0,nn)! Multiples of r.l.v. that bound cutoff gmax0  (r.l.v. reciplocal lattice vector)
  implicit none
  intent(in) ::   q0,q1,q2,qp,mshlst,gmax0
  intent(out)::                            nn
  ! ----------------------------------------------------------------------
  !i Inputs
  !i   q0    :r.l.v. for which bounds nmin and nmax are to be computed
  !i   q1    :first  r.l.v. different from q0
  !i   q2    :second r.l.v. different from q0
  !i   qp    :k-point added to G vectors: sphere is centered G=qp.
  !i   mshlst:An ordered list of integers used to restrict the assignment
  !i         :of nn to one of an allowed list of points, should a value
  !i         :of nn be assigned (see description for nn below).
  !i         :If first entry is nonzero, mshlst(1..) = list of allowed
  !i         :values nn may take.  mshlst(0) is the size of mshlst.
  !i         :Certain fourier transform routines have restrictions
  !i         :on the allowed mesh sizes; this constraint is designed
  !i         :to handle that restriction.
  !i  gmax0  :cutoff G
  !o Outputs
  !o     nn :   meshsize
  !r Remarks
  !r   q0,q1,q2,qp and G are all dimensionless (units of 2*pi/a)
  ! ----------------------------------------------------------------------
  integer :: nmin,nmax,nn,mshlst(0:*), indx,i
  double precision :: q0(3),q1(3),q2(3),qp(3),gmax0,qperp(3),ddot,qqperp
  ! ... qperp = q1 x q2 / |q1 x q2| ; qqperp = q . qperp
  qperp=[q1(2)*q2(3) - q1(3)*q2(2),&
       q1(3)*q2(1) - q1(1)*q2(3),&
       q1(1)*q2(2) - q1(2)*q2(1)]
  qperp  = qperp/sum(qperp**2)**.5 
  qqperp = sum(q0*qperp)
  nmax =  gmax0/abs(Qqperp) - sum(qp*qperp)/Qqperp + 1
  nmin = -gmax0/abs(Qqperp) - sum(qp*qperp)/Qqperp - 1
  nn = 2*max(iabs(nmin),nmax)+1
  if (mshlst(0) /= 0) then
     indx = 1
     !call hunti(mshlst(1),mshlst,nn,0,indx)
     indx= findloc([(mshlst(i)<nn.and.nn<=mshlst(i+1),i=1,mshlst(0)-1)],value=.true.,dim=1)
     nn = mshlst(min(indx+1,mshlst(0)))
  endif
end subroutine gvlstn
subroutine ppfac(fmax,job,fac,nfac)   !- Find all products of prime factors within some maximum
  !i Inputs
  !i   npfac :number of prime factors
  !i   pfac  :list of prime factors
  !i   fmax  :maximum value of product
  !i   job   :0 list is returned unsorted
  !i         :1 list is returned sorted
  !o Outputs
  !o   fac   :all products of pfac <= fmax
  !o         :List is returned sorted.
  !o         :fac must be dimensioned at least nfac
  !o         :or 3*nfac if fac is returned sorted.
  !o   nfac  :length of fac
  ! ----------------------------------------------------------------------
  implicit none
  integer,parameter:: npfac=5,pfac(1:npfac)=[2,3,5,7,11]
  integer :: fac(*),fmax,job,nfac
  integer :: i,m,k,fack,nfaci
  nfac = 0
  nfaci = 0
  do  i = 1, npfac
     nfac = nfac+1
     nfaci = nfac
     fac(nfac) = pfac(i)
     do  k  = 1, nfaci
        fack = fac(k)
        do  m = 1, fmax
           fack = fack*pfac(i)
           if (fack <= fmax) then
              nfac = nfac+1
              fac(nfac) = fack
           else
              goto 10
           endif
        enddo
10      continue
     enddo
  enddo
  if (job /= 0) call ivheap(1,nfac,fac,fac(1+2*nfac),0)
end subroutine ppfac
! subroutine hunti(xa,n,x,iprm,low)  !- Brackets a value within an ordered array of integer points
!   ! ----------------------------------------------------------------
!   !i Inputs
!   !i   xa  :array of points
!   !i   n   :size of array
!   !i   x   :value to bracket
!   !i   iprm:permutation table by which array xa is ordered
!   !i        iprm(1) <= 0 => assumes iprm(i) = i; iprm not referenced
!   !i   low : initial guess for output low
!   !o Outputs
!   !o   low : xa(low) < x <= xa(low+1)
!   !u Updates
!   !u   13 Sep 01 Handle case n=0
!   ! ----------------------------------------------------------------
!   !     implicit none
!   ! Passed parameters
!   integer :: n,low,iprm(n)
!   integer :: xa(n),x
!   ! Local variables
!   integer :: inc,jhi,jm
!   logical :: ascnd,liprm
!   integer :: xn

!   if (n == 0) then
!      low = 0
!      return
!   endif

!   liprm = iprm(1) .gt. 0
!   if (liprm) then
!      ascnd = xa(iprm(n)) .gt. xa(iprm(1))
!   else
!      ascnd = xa(n) .gt. xa(1)
!   endif
!   if (low <= 0 .OR. low > n) then
!      low = 0
!      jhi = n+1
!      goto 3
!   endif
!   inc = 1
!   if (liprm) then
!      xn = xa(iprm(low))
!   else
!      xn = xa(low)
!   endif
!   if (x >= xn .eqv. ascnd) then
! 1    jhi = low+inc
!      if (jhi > n) then
!         jhi = n+1
!      else
!         if (liprm) then
!            xn = xa(iprm(jhi))
!         else
!            xn = xa(jhi)
!         endif
!         if (x >= xn .eqv. ascnd) then
!            low = jhi
!            inc = inc+inc
!            goto 1
!         endif
!      endif
!   else
!      jhi = low
! 2    low = jhi-inc
!      if (low < 1) then
!         low = 0
!      else
!         if (liprm) then
!            xn = xa(iprm(low))
!         else
!            xn = xa(low)
!         endif
!         if (x < xn .eqv. ascnd) then
!            jhi = low
!            inc = inc+inc
!            goto 2
!         endif
!      endif
!   endif
! 3 if (jhi-low == 1) then
!      !   ... Find the first of values equal to x
! 4    continue
!      if (low > 1) then
!         if (liprm) then
!            xn = xa(iprm(low-1))
!         else
!            xn = xa(low-1)
!         endif
!         if (xn == x) then
!            low = low-1
!            goto 4
!         endif
!      endif
!      !     ... if xa(low) = x, decrement low
!      if (low >= 1) then
!         if (liprm) then
!            xn = xa(iprm(low))
!         else
!            xn = xa(low)
!         endif
!         if (xn == x) low = low-1
!      endif
!      return
!   endif
!   jm = (jhi+low)/2
!   if (liprm) then
!      xn = xa(iprm(jm))
!   else
!      xn = xa(jm)
!   endif
!   if (x > xn .eqv. ascnd) then
!      low = jm
!   else
!      jhi = jm
!   endif
!   goto 3
! end subroutine hunti


