subroutine gridls(x,y,sigmay,npts,nterms,mode,a,deltaa,
     & sigmaa,yfit,chisqr,ra,dec,maxp)
c
c     This subroutine is a slightly modified version of GRIDLS in
c     "Data Reduction and Error Analysis for the Physical Sciences" by
c     P.R. Bevington (1969), p.212
c
      implicit undefined(a-z)
      integer*4 maxp,nterms
      real*8 x(maxp),y(maxp),sigmay(maxp),yfit(maxp),ra(maxp),dec(maxp)
      real*8 a(nterms),deltaa(nterms),sigmaa(nterms)
      real*8 chisqr,chisq1,chisq2,chisq3
      real*8 delta,save
      real*4 fn
      integer*4 nfree,npts,mode,j
      character text*80
      real*8 fchisq
      nfree=npts-nterms
      chisqr=0.
      if(nfree.le.0)goto 100
      do 90 j=1,nterms
c
c     Evaluate chi square at first two search points
c
d     write(6,*)' GRIDLS:  going to FIXELL first time'
      call fixell(x,y,npts,a,yfit,ra,dec,maxp)
d     write(6,*)' GRIDLS:  going to FCHISQ first time'
      chisq1=fchisq(y,sigmay,npts,nfree,mode,yfit,maxp)
d     write(6,'(a,f20.10)')' GRIDLS:  chisq1=',chisq1
      fn=0.
      delta=deltaa(j)
   41 a(j)=a(j)+delta
d     write(6,*)' GRIDLS:  going to FIXELL second time'
      call fixell(x,y,npts,a,yfit,ra,dec,maxp)
d     write(6,*)' GRIDLS:  going to FCHISQ second time'
      chisq2=fchisq(y,sigmay,npts,nfree,mode,yfit,maxp)
d     write(6,'(a,f20.10)')' GRIDLS:  chisq2=',chisq2
      if(chisq1-chisq2.eq.0)goto 41
      if(chisq1-chisq2.gt.0)goto 61
c
c     reverse direction of search if chi square increases
c
      delta=-delta
      a(j)=a(j)+delta
      save=chisq1
      chisq1=chisq2
      chisq2=save
c
c     increment a(j) until chi square increases
c
   61 fn=fn+1.
      if(fn.gt.50.)then
       write(text,1)j
    1  format('Model parameter',i2,' did not converge',
     & ' in 50 iterations!!')
       call putout(text)
       return
      endif
      a(j)=a(j)+delta
d     write(6,*)' GRIDLS:  going to FIXELL third time'
      call fixell(x,y,npts,a,yfit,ra,dec,maxp)
d     write(6,*)' GRIDLS:  going to FCHISQ third time'
      chisq3=fchisq(y,sigmay,npts,nfree,mode,yfit,maxp)
d     write(6,'(a,f20.10)')' GRIDLS:  chisq3=',chisq3
      if(chisq3-chisq2.ge.0)goto 81
      chisq1=chisq2
      chisq2=chisq3
      goto 61
c
c     find minimum of parabola defined by last three points
c
   81 delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5)
      a(j)=a(j)-delta
      sigmaa(j)=deltaa(j)
     &   *sqrt(2./(dble(nfree)*(chisq3-2.*chisq2+chisq1)))
      deltaa(j)=deltaa(j)*fn/3.
   90 continue
c
c     evaluate fit and chi square for final parameters
c
d     write(6,*)' GRIDLS:  going to FIXELL fourth time'
      call fixell(x,y,npts,a,yfit,ra,dec,maxp)
d     write(6,*)' GRIDLS:  going to FCHISQ fourth time'
      chisqr=fchisq(y,sigmay,npts,nfree,mode,yfit,maxp)
d     write(6,'(a,f20.10)')' GRIDLS:  chisqr=',chisqr
  100 return
      end