c     subroutine ascent(xp,yp,zp,origin,stepc,val,niter,ir)
      subroutine ascent(xp,yp,zp,val,niter,ir)
      use grid_module
      use basins_module
      implicit none
      double precision xp,yp,zp,val,xold,yold,zold,valold,stepc
      double precision rho,drho(3),d2rho(6),elfd,gelf(3),d2elf(6)
      double precision grad_norm,delta,df(3),f,xv,yv,zv,dist
      double precision deltaold,grad_old
      double precision gold(3),tilt,rand
      integer niter,level,icell(8),iacc
      integer i,j,ixa,iya,iza,indy,npyz,ir,isave
      logical found,func_type,test_cell
      common/ftype/func_type
      common/accuracy/iacc
      xold=xp
      yold=yp
      zold=zp
      xv=xp
      yv=yp
      zv=zp
      isave=iacc
      npyz=np(2)*np(3)
      delta=0.0d0
      found=.true.
      niter=0
      ir=0
      level=1
      gold(1)=0.1d0
      gold(2)=0.0d0
      gold(3)=0.0d0
      xp=xp+0.005d0*(1.0d0-rand)*step(1)
      yp=yp+0.005d0*(1.0d0-rand)*step(2)
      zp=zp+0.005d0*(1.0d0-rand)*step(3)
      stepc=0.0625d0*dsqrt(step(1)**2+step(2)**2+step(3)**2)
      if(iacc.eq.0) then
         call elf_eval(rho,drho,d2rho,elfd,gelf,d2elf,xp,yp,zp,level)
         if(func_type) then
            df(1)=gelf(1)
            df(2)=gelf(2)
            df(3)=gelf(3)
            f=elfd
         else
            df(1)=drho(1)
            df(2)=drho(2)
            df(3)=drho(3)
            f=rho
         endif
      else
         if(func_type) then 
            call elf_dif(f,df,xp,yp,zp)
         else
            call rho_dif(f,df,xp,yp,zp)
         endif
      endif
      grad_old=gold(1)**2+gold(2)**2+gold(3)**2
      grad_old=dsqrt(grad_old)
      grad_norm=grad_old
      gold(1)=df(1)/grad_old
      gold(2)=df(2)/grad_old
      gold(3)=df(3)/grad_old
      valold=f
      do while(found)
         niter=niter+1
         grad_old=grad_norm
         grad_norm=df(1)**2+df(2)**2+df(3)**2
         grad_norm=dsqrt(grad_norm)
         tilt=df(1)*gold(1)+df(2)*gold(2)+df(3)*gold(3)
         tilt=tilt/grad_norm
         xp=xold+stepc*df(1)/grad_norm
         yp=yold+stepc*df(2)/grad_norm
         zp=zold+stepc*df(3)/grad_norm
         if(step(1).gt.0.0d0) then
            ixa=int((xp-origin(1))/step(1))
         else
            ixa=1
         endif
         if(step(2).gt.0.0d0) then
            iya=int((yp-origin(2))/step(2))
         else
            iya=1
         endif
         if(step(3).gt.0.0d0) then
            iza=int((zp-origin(3))/step(3))
         else
            iza=1
         endif
         if((ixa.ge.0).and.(ixa.lt.np(1))) then
         if((iya.ge.0).and.(iya.lt.np(2))) then
         if((iza.ge.0).and.(iza.lt.np(3))) then
         indy=1+iza+np(3)*(iya+np(2)*ixa)
         icell(1)=attract_code(indy)
         icell(2)=attract_code(indy+1)
         icell(3)=attract_code(indy+np(3))
         icell(4)=attract_code(indy+np(3)+1)
         icell(5)=attract_code(indy+npyz)
         icell(6)=attract_code(indy+npyz+1)
         icell(7)=attract_code(indy+npyz+np(3))
         icell(8)=attract_code(indy+npyz+np(3)+1)
         test_cell=.true.
         do i=2,8
            if(icell(i).ne.icell(1)) then
               test_cell=.false.
            endif
         enddo
         endif
         endif
         endif
         if(ir.eq.0) test_cell=.false.
         if(test_cell) then
            ir=icell(1)
            return
         endif
         gold(1)=df(1)/grad_norm
         gold(2)=df(2)/grad_norm
         gold(3)=df(3)/grad_norm
         do j=1,attract_label
            dist=(xp-x(j))**2+(yp-y(j))**2+(zp-z(j))**2
            dist=dsqrt(dist)
            if(dist.le.0.25d0) then
               xp=x(j)
               yp=y(j)
               zp=z(j)
               val=val_elf(j)
               ir=j
               return
            endif
         enddo
         if((f.lt.valold).or.(tilt.le.0.5d0)) then
           stepc=stepc*0.3819660
           xp=xold
           yp=yold
           zp=zold
         else
           deltaold=delta
           delta=f-valold
           xold=xp
           yold=yp
           zold=zp
           if(grad_norm.le.grad_old) then
              stepc=stepc*0.9819660
           else if((stepc.lt.0.2d0).and.(tilt.ge.0.8d0)) then
              stepc=stepc/0.3819660
           endif
         endif
         dist=(xp-xv)**2+(yp-yv)**2+(zp-zv)**2
         if(dist.gt.1.0d2) then
            found=.false.
            valold=0.0d0
         endif
         if(grad_norm.lt.1.0d-6) then
            found=.false.
         endif
         if(stepc.lt.1.0d-4) then
            found=.false.
         endif
         if(niter.gt.50) found=.false.
         valold=f
         if(iacc.eq.0) then
            call elf_eval(rho,drho,d2rho,elfd,gelf,d2elf,xp,yp,zp,level)
            if(func_type) then
               f=elfd
               df(1)=gelf(1)
               df(2)=gelf(2)
               df(3)=gelf(3)
            else
               f=rho
               df(1)=drho(1)
               df(2)=drho(2)
               df(3)=drho(3)
            endif
         else
            if(func_type) then
               call elf_dif(f,df,xp,yp,zp)
            else
               call rho_dif(f,df,xp,yp,zp)
            endif
         endif
      enddo
      ir=0
      return
      end
