C------ CONFUN & MAKEOPT are for use by NPSOL, otherwise not used 
        SUBROUTINE CONFUN(MODE,NCNLN,N,NROWJ,NEEDC,X,C,CJAC,NSTATE)
        IMPLICIT double precision (A-H,O-Z)
        INTEGER*4 MODE,NCNLN,N,NROWJ
        INTEGER*4 NEEDC(*)
        REAL*8 X(N),C(*),CJAC(NROWJ,*)
        SAVE
        RETURN
        END

        subroutine makeopt(bigbnd,itmax1,itmax2,ftol1,ftol2,msglvl,
     *                     mverify)
        implicit double precision(a-h,o-z)

      OPEN(7,FILE='options1.dat',status='unknown')
      rewind(7)
      WRITE(7,7000)
      write(7,7010)
      WRITE(7,7001) BIGBND
      WRITE(7,7002) ITMAX1
      WRITE(7,7004) FTOL1
      WRITE(7,7005) MSGLVL
      WRITE(7,7006) MVERIFY
      WRITE(7,7900)
      close(7)


      OPEN(7,FILE='options2.dat',status='unknown')
      rewind(7)
      WRITE(7,7000)
      write(7,7010)
      WRITE(7,7001) BIGBND
      WRITE(7,7002) ITMAX2
      WRITE(7,7004) FTOL2
      WRITE(7,7005) MSGLVL
      WRITE(7,7006) MVERIFY
      WRITE(7,7900)
      close(7)

      ITMAX3=100
      OPEN(7,FILE='options3.dat',status='unknown')
      rewind(7)
      WRITE(7,7000)
      write(7,7010)
      WRITE(7,7001) BIGBND
      WRITE(7,7002) ITMAX3
      WRITE(7,7004) FTOL2
      WRITE(7,7005) MSGLVL
      WRITE(7,7006) MVERIFY
      WRITE(7,7900)
      close(7)

7000  FORMAT('BEGIN',74X,' ')
7010  format('   NOLIST                    ',48x,' ')
7001  FORMAT('   INFINITE BOUND            ',D15.3,35X,' ')
7002  FORMAT('   MAJOR ITERATION LIMIT     ',I15  ,35X,' ')
7004  FORMAT('   OPTIMALITY TOLERANCE      ',D15.3,35X,' ')
7005  FORMAT('   MAJOR PRINT LEVEL         ',I15  ,35X,' ')
7006  FORMAT('   VERIFY LEVEL              ',I15  ,35X,' ')
7008  FORMAT('   LINESEARCH TOLERANCE      ',D15.3,35X,' ')
7900  FORMAT('END  ',74X,' ')
      return
      end


        subroutine sdev(x,n1,n,sd,ave)
c
c    Purpose: calculates standard deviation of values n1,n2,....n
c             in array x of length n
c
        implicit double precision(a-h,o-z)
        dimension x(n)
        xn=float(n-n1+1)
        ave=0.
        do 10 i=n1,n
            ave=ave+x(i)
10      continue
        ave=ave/xn
        var=0.
        do 20 i=n1,n
            var=var+(x(i)-ave)**2
20      continue
        var=var/(xn-1.)
        sd=sqrt(var)
        return
        end

        integer*4 function igcdiv(i1,i2)
        implicit integer*4(i-n)
        igcd=0
        imin=min(i1,i2)
        do 10 i=1,imin
            ir1=mod(i1,i)
            ir2=mod(i2,i)
            if (ir1.eq.0.and.ir2.eq.0) igcd=i
10      continue
        igcdiv=igcd
        return
        end

       subroutine getx(fname,xdat,nxmax,nx,ifreq,ltr)
       implicit double precision(a-h,o-z)
       character*30 fname
       dimension xdat(nxmax)
       open(unit=11,file=fname,status='old')
       rewind(11)
       iskip=ifreq-1
       m=0
       do 50 i=1,nx
           read(11,*,end=60) xdat(i)
           m=m+1
           if (ifreq.gt.1) then
                do 40 j=1,iskip
                    read(11,*,end=60) xjunk
40              continue
           endif
           if (ltr.gt.0) xdat(i)=dlog(xdat(i))
50      continue
60     continue
       close(11)
       if (m.lt.nx) nx=m
       return
       end

        subroutine recontp(xt,xdat,xlag,nxmax,jdmax,nx,ldelay,
     *     jd,jt1,m,jtp)
        implicit double precision(a-h,o-z)
        dimension xt(nxmax),xdat(nxmax),xlag(nxmax,jdmax)
        m=0
        do 10 i=1,nx-jt1
           m=m+1
           ipred=i+jt1
           xt(i)=xdat(ipred)
	   if (jd.gt.0) then	
             do 5 j=1,jd
               iback=jtp+(j-1)*ldelay
               xlag(i,j)=xdat(ipred-iback)
5            continue
	   endif
10      continue
        return
        end


      double precision function uran(iseed)
      integer  iy
c    ==> the uran code was grabbed from netlib at ORNL.
c
c      urand is a uniform random number generator based  on  theory  and
c  suggestions  given  in  d.e. knuth (1969),  vol  2.   the integer  iy
c  should be initialized to an arbitrary integer prior to the first call
c  to urand.  the calling program should  not  alter the  value  of  iy
c  between  subsequent calls to urand. values of urand will be returned
c  in the interval (0,1).
c
c    ==> uran is modified for LENNS so that if iseed.le.0, the generator is
c        reseeded, otherwise it iterates from the current value of iy
c        (stored internally), so the calling program can use any
c        positive value of iseed after the generator has been seeded.
c
      integer  ia,ic,itwo,m2,m,mic
      double precision  halfm,uran,s
      double precision  datan,dsqrt
      data m2/0/,itwo/2/
      if(iseed.le.0) iy=iseed
      if (m2 .ne. 0) go to 20
c
c  if first entry, compute machine integer word length
c
      m = 1
   10 m2 = m
      m = itwo*m2
      if (m .gt. m2) go to 10
      halfm = m2
c
c  compute multiplier and increment for linear congruential method
c
      ia = 8*idint(halfm*datan(1.d0)/8.d0) + 5
      ic = 2*idint(halfm*(0.5d0-dsqrt(3.d0)/6.d0)) + 1
      mic = (m2 - ic) + m2
c
c  s is the scale factor for converting to floating point
c
      s = 0.5/halfm
c
c  compute next random number
c
   20 iy = iy*ia
c
c  the following statement is for computers which do not allow
c  integer overflow on addition
c
      if (iy .gt. mic) iy = (iy - m2) - m2
c
      iy = iy + ic
c
c  the following statement is for computers where the
c  word length for addition is greater than for multiplication
c
      if (iy/2 .gt. m2) iy = (iy - m2) - m2
c
c  the following statement is for computers where integer
c  overflow affects the sign bit
c
      if (iy .lt. 0) iy = (iy + m2) + m2
      uran = float(iy)*s
      return
      end