subroutine gridls(x,y,sigmay,npts,nterms,mode,a,deltaa,
     & sigmaa,yfit,chisqr)
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
      real*8 x(*),y(*),sigmay(*),yfit(*)
      real*8 a(*),deltaa(*),sigmaa(*)
      real*8 chisqr,chisq1,chisq2,chisq3,delta,save

      modef=0 ! Gaussian model 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
      if(mode.eq.0)then
       call rtfuncs(npts,x,yfit,a)
      elseif(mode.eq.1)then
       call fixell(x,y,npts,a,yfit)
      else
       write(6,*)' GRIDLS: mode not allowed:',mode
       stop
      endif
      call fchisq(y,sigmay,npts,nfree,yfit,chisq1,modef)
      fn=0.
      delta=deltaa(j)
   41 a(j)=a(j)+delta
      if(mode.eq.0)then
       call rtfuncs(npts,x,yfit,a)
      elseif(mode.eq.1)then
       call fixell(x,y,npts,a,yfit)
      endif
      call fchisq(y,sigmay,npts,nfree,yfit,chisq2,modef)
      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.100)then
       write(6,*)'GRIDLS: Model parameter',j,' did not converge',
     & ' in 100 iterations!!'
       write(6,*)'a(',j,')=',a(j),' delta=',delta
       stop
      endif
      a(j)=a(j)+delta
      if(mode.eq.0)then
       call rtfuncs(npts,x,yfit,a)
      elseif(mode.eq.1)then
       call fixell(x,y,npts,a,yfit)
      endif
      call fchisq(y,sigmay,npts,nfree,yfit,chisq3,modef)
      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.D0/(1.D0+(chisq1-chisq2)/(chisq3-chisq2))+0.5D0)
      a(j)=a(j)-delta
      sigmaa(j)=deltaa(j)
     &   *sqrt(2.D0/(dble(nfree)*(chisq3-2.D0*chisq2+chisq1)))
      deltaa(j)=deltaa(j)*dble(fn/3.)
   90 continue
c
c     evaluate fit and chi square for final parameters
c
      if(mode.eq.0)then
       call rtfuncs(npts,x,yfit,a)
      elseif(mode.eq.1)then
       call fixell(x,y,npts,a,yfit)
      endif
      call fchisq(y,sigmay,npts,nfree,yfit,chisqr,modef)
  100 return
      end