      program lenns
      implicit double precision(a-h,o-z)
      parameter(kmax=10,jdmax=10,nxmax=525,lmax=65)
      parameter(npmax=1+kmax*(jdmax+2),twopi=6.2831853)
      dimension xt(nxmax),xdat(nxmax),xlag(nxmax,jdmax),
     *    theta(npmax),svexp(lmax),qrexp(lmax),
     *    der(lmax,nxmax),tsav(20,npmax),rsav(20),xraw(nxmax)
      dimension h0(npmax,npmax)
      character*35 fname,outname
      integer tim1,tim2,time
      common /size/ k,jd,nx,m,jforce
      common /xdata/ xt,xlag
      external objfun

      tim1=time()
c----------- Initialize random number generator
      iseed=time()
      rn=uran(-abs(mod(iseed,15000)))

c----------- Set initial Jacobian to identity matrix
      call setidmat(h0,npmax)

c----------- Read in parameters and minimization routine options
      open(11,file='lenns.par',status='old')
      read(11,10) fname
      read(11,10) outname
10    format(a35)
      read(11,*) nx,jt1,ltr,ifreq,jd1,jd2,k1,k2,kcrit,l1,l2,
     *     jtpsav,alf,jforce,jper,cgcv,scale,ftol1,ftol2,
     *     itmax1,itmax2,maxstep1,maxstep2
      close(11)
      iprint=itmax1+itmax2

c----------- Check if jt1 is OK, and fix if necessary
      if (jtpsav.lt.0) then
          jt1min=jd2*l2
      else
          jt1min=(jd2-1)*l2+jtpsav
      endif
      if (jt1.lt.jt1min) jt1=jt1min

c------------ read data series into xdat; set nx=total series length
      call getx(fname,xdat,nxmax,nx,ifreq,ltr)

c---------- Set up the output file
      call makeout(fname,outname,nx,jt1,ltr,ifreq,k1,k2,cgcv,
     *     kcrit,jtpsav,alf,scale,jforce,jper,ftol1,ftol2,
     *     itmax1,itmax2,maxstep1,maxstep2)

c------------ Standardize values to predict: SD=1, mean=0 if alf=0
      call sdev(xdat,jt1+1,nx,sd,xbar)
      if (alf.eq.0.) then
          do 150 j=1,nx
              xdat(j)=xdat(j)-xbar
150       continue
      endif
      call dscal(nx,1./sd,xdat,1)

c------------------ START MAIN LOOP on time-delay and #lags
      do 1010 ldelay=l1,l2
          do 1005 jd=jd1,jd2

              jtp=jtpsav
              if (jtpsav.le.0) jtp=ldelay

              call recontp(xraw,xdat,xlag,nxmax,jdmax,nx,ldelay,jd,jt1,
     *            m,jtp)
c --------------- m is the # of state vectors (# rows in xlag)
              do 200 j=1,m
                  xt(j)=xraw(j)-alf*xlag(j,1)
200           continue

c----------------- Add sin & cos terms if fitting forced model
              if (jforce.gt.0) then
                  xper=float(jper)
                  do 210 j=1,m
                      ajper=float(mod(j,jper))/xper
                      xlag(j,jd+1)=cos(twopi*ajper)
                      xlag(j,jd+2)=sin(twopi*ajper)
210               continue
              endif

c---------------- Initialize markers for best # units 'so far'
              kbic=k1
              bicopt=10000.
              kgcv=k1
              gcvopt=10000.
              mbic=0
              mgcv=0
c------------------ START LOOP on # units
              do 1004 k=k1,k2

c ---------------- decide if it's time to quit increasing # units
                  if (k.gt.(2+kbic)) mbic=mbic+1
                  if (k.gt.(2+kgcv)) mgcv=mgcv+1
                  if (kcrit.eq.0) kquit=mbic
                  if (kcrit.eq.1) kquit=mgcv
                  if (kcrit.eq.2) kquit=mbic*mgcv
                  if (kquit.gt.0) goto 1004
                  npar=1+k*(jd+2+2*jforce)
                  if (2*npar.gt.m) goto 1004
                  if (npar.gt.npmax) goto 1004

                  do 15 i=1,20
                      rsav(i)=10000.
15                continue

c------------------ START LOOP on repeated fits
                  do 99 irep=1,251

c----------------- initialize parameter choices
                      eps=scale*10**( .01*(dfloat(irep)-126.) )
                      call tinit(theta,npar,eps,100)

c------------------ do a high-tolerance fit
                      call bfgsfm(objfun,npar,npmax,theta,h0,ftol1,
     *                    itmax1,maxstep1,fnew,iter,iprint,inform)

c                      call npsol(npar,nclin,ncnln,nrowa,nrowj,nrowr,a,bl,bu,
c     *                    confun,objfun,inform,iter,istate,c,cjac,clamda,
c     *                    fnew,g,r,theta,iw,leniw,w,lenw)


                      xm=float(m)
                      bic=1.419+dlog(fnew)+0.5*float(npar)*dlog(xm)/xm
                      imod=mod(irep-1,10)
                      if(imod.eq.0) then
                          open (13,file='lenns.tmp',status='unknown',
     *                        access='append')
                       write(13,97) ldelay,jd,k,eps,fnew,bic,iter,inform
                          close(13)
                      endif
                      imax=idamax(20,rsav,1)
                      rmax=rsav(imax)
                      if (fnew.lt.rmax) then
                          rsav(imax)=fnew
                          do 85 jpar=1,npar
                              tsav(imax,jpar)=theta(jpar)
85                        continue
                      endif
99                continue
c---------------------- END LOOP on repeated fits

                  do 999 i=1,20
                      do 997 jpar=1,npar
                          theta(jpar)=tsav(i,jpar)
997                   continue

                      call bfgsfm(objfun,npar,npmax,theta,h0,ftol2,
     *                    itmax2,maxstep2,fnew,iter,iprint,inform)

c                      call npsol(npar,nclin,ncnln,nrowa,nrowj,nrowr,a,bl,bu,
c     *                   confun,objfun,inform,iter,istate,c,cjac,clamda,
c     *                   fnew,g,r,theta,iw,leniw,w,lenw)
              (7,g,r,theta,iw,leniw,w,lenw)
    