	program rcs
c
c	reads an raw ring width increment file and prints
c	out a list of all records held within
c
c	pjk_trl may 12, 1992, modified 11/2000
c
c
	integer mxdat,iu,ou,mxfl,mxrec
	parameter(mxdat=500,mxrec=50,mxfl=40,
     *           iu=13,ou=15)
       parameter(flag=-999.0)
	integer fnum,sumrrw,eof,p,j,k,m,l,i,itype
	integer iyr(mxrec),lyr(mxrec),age(mxrec),nc(mxrec)
	real*4    rrw(mxdat),sum,mrw
	character*8 id(mxrec)
	character*40 infl(mxfl),oufl
	character*1 x
	character*2 opt
	
	real*4 xrc(mxdat),f(mxdat),wk2(mxdat),wk3(mxdat),eq(7)
	
	real*4 x1(mxdat,mxrec)
	integer idate(mxrec,3),pith(mxrec),jdt,ndt,n1
	logical ok1,ok2
	logical poff
c
c	begin program
c
	ok1=.false.
c
	write(*,'(/13x,''Program MakeRec_v10 reports the contents of'')')
	write(*,'(13x,''of decade format Ring-Width data files and/or '')')
	write(*,'(13x ''ARSTAN Chronology files.'')')
	write(*,'(/13x,''Current dimentions are:'')')
	write(*,'(13x,''  Max: records in a single file = '',i8)')mxrec
	write(*,'(13x,''  Max: record length            = '',i8)')mxdat
	write(*,'(13x,''  Max: files                    = '',i8,/)')mxfl
c	write(*,'(/13x,''Read in file names from a .lst file?'')')
c	write(*,'(13x,''  [This is a simple text file containing one '')')
c	write(*,'(13x,''  column of data file names.]'',/)')
c	write(*,'(13x,''Answer [T]rue or [F]alse :  '',$)')
c	read(*,*)ok1
	ok1=.false.
c
	if(.not.ok1)then
	  write(*,' (/13x,''Enter file name :  '',$)')
	  read(*,'(a)')infl(mxfl)
	  ok2=.true.
	  inquire(file=infl(mxfl),exist=ok2)
	  if(.not.ok2)then
	     write(*,*)'Warning'
	     write(*,'(/1x,''This file '',a,'' DOES NOT exist'',/)')infl(j)
	     write(ou,*)'Warning'
	     write(ou,'(/1x,''This file '',a,'' DOES NOT exist'',/)')infl(j)
	     write(*,'(1x,''Program Terminating'')')
	     pause 'Hit RETURN to Quit'
	     stop
	  endif
	  do j=1,mxfl-1
	     infl(j)=infl(mxfl)
	  end do
	  fnum=1
	  write(*,*)
	  open(unit=iu,file=infl(mxfl),status='old',readonly)
	else
	  write(*,'(/13x,"List- file name ==> ",$)')
	  read(*,'(a)')infl(mxfl)
	  open(unit=iu,file=infl(mxfl),status='old',readonly)
	  write(*,'(/13x,''Openning files:'')')
	  fnum=0             
	  do p=1,mxfl-1
5	     read(iu,'(a)',end=7)infl(p)
	     write(*,'(13x,i4,'' '',a)')p,infl(p)
	     fnum=fnum+1
	  end do
7	  close(iu)
	endif	
c
	do p=1,40
         read(infl(mxfl)(p:p),'(a1)')x
	  if(x.eq.' ')goto 15
	end do
15	p=p-1
	oufl=infl(mxfl)(1:p)//'_rec'
	open(unit=ou,file=oufl,status='unknown')
	write(ou,'(/1x,''Input file for this run : '',a)')infl(mxfl)
	write(*,'(/13x,''Input file for this run : '',a)')infl(mxfl)
	write(*,'(/13x,''Measurements or Chronologies?'')')
	write(*,'(13x,''[1] rrw  [2] crn  : '',$)')
	read(*,'(i1)')itype
	write(*,*)
	if(itype.eq.1)opt='mr'
	if(itype.eq.2)opt='ir'
c
	sumrrw=0
	p=1
	do j = 1,fnum      ! loop to read all files in list file
	  k=1
	  write(ou,'(1x'' ___ '',a,t42,'' ___ '')')infl(j)
	  write(*,'(1x,'' ___ '',a,t42,'' ___ '')')infl(j)
	  if(ok1)then
	     ok2=.true.
	     inquire(file=infl(j),exist=ok2)
	     if(.not.ok2)then
	        write(*,*)'Warning'
	        write(*,'(/1x,''This file '',a,'' DOES NOT exist'',/)')infl(mxfl)
	        write(ou,*)'Warning'
	        write(ou,'(/1x,''This file '',a,'' DOES NOT exist'',/)')infl(mxfl)
	        goto 30
	     else
	        open(unit=iu,file=infl(j),status='old',readonly)
	     endif
	  endif
	  if(itype.eq.1)then
	    write(ou,'(1x,''        ID        IYR     LYR     N      m_rrw'')')
	  else
	    write(ou,'(1x,''        ID        IYR     LYR     N      m_indx'')')
	  endif
20       call trrw(id(p),iyr(p),age(p),rrw,nc,eof,iu,6,opt)
	     sum=0.0
	     do m=1,age(p)
	        sum=sum+rrw(m)
	        x1(m,p)=rrw(m)
	     enddo
	     mrw=sum/age(j)
c
	     if(eof.gt.0)then
	        lyr(p)=iyr(p)+age(p)-1
	        write(ou,'(1x,i8,1x,a,3(i8),f8.3)')k,id(p),iyr(p),lyr(p),age(p),mrw
	        sumrrw=sumrrw+age(p)
	        pith(p)=p; idate(p,1)=iyr(p); idate(p,2)=lyr(p);idate(p,3)=age(p)
	        p=p+1
	        k=k+1
	        goto 20
	     else
	        close(iu)
	        write(*,'(/2x,'' Series Read from file= '',i8)')k-1
	        write(*,'(2x,'' Total Series Read= '',i8)')p-1
	        eof=0
	     endif
30	     continue
	end do
c
	nc(1)=p-1
c
	if(nc(1).le.mxrec)call comint(mxrec,iyr,id,age,nc(1),ou)	        
	close(ou)
	if(.not.ok1)close(iu)
	
40          write(*,'(/1x,''you can change the current curve '',
     *      ''fit with any one of the following options:''/)')
          write(*,'(t3,''  -5:  friedman super smoother'',
     *      ''        lower alpha (0-9) gives more local fit'')')
          write(*,'(t3,''  -4:  robust lowess smoother'',
     *      ''         integer percent of series length'')')
          write(*,'(t3,''  -3:  median smoothing'',
     *      ''               must have odd window length'')')
          write(*,'(t3,''   1:  neg expon curve, '',
     *      ''no = opt 4    f(i) = a*exp(-b*t(i)) + d    (d>0)'')')
          write(*,'(t3,''   2:  neg expon curve, '',
     *      ''no = opt 5    f(i) = a*exp(-b*t(i)) + d    (d>0)'')')
          write(*,'(t3,''   3:  neg expon curve  '',
     *      ''(general)     f(i) = a*exp(-b*t(i+p)) + d (any d)'')')
          write(*,'(t3,''   4:  linear regression '',
     *      ''(any slope)  f(i) = +/-c*t(i) + d'')')
          write(*,'(t3,''   5:  linear regression '',
     *      ''(neg slope)  f(i) = -c*t(i) + d   '')')
          write(*,'(t3,''   6:  horizontal line through '',
     *      ''mean   f(i) = mean(y(i)) = d'')')
          write(*,'(t3,''   7:  hugershoff growth '',
     *      ''function     f(i) = a*t(i+p)**b * exp(-c*t(i+p)) + d'')')
          write(*,'(t3,''   8:  general exponential (b = 1)    '',
     *      ''f(i) = a*t(i+p)**b * exp(-c*t(i+p)) + d'')')
          write(*,'(t3,''  >9: cubic smoothing '',
     *      ''spline         fixed 50 pct variance cutoff'')')
          write(*,'(t3,'' <-9: cubic smoothing '',
     *      ''spline         pct n 50 pct variance cutoff'')')
               write(*,'(/1x,''enter new option (h for help)'',/
     *      '' ==> '',$)')
       read(*,*)jdt
	write(*,*)
	
c make a dummy pith offset record
	do p=1,nc(1)
	  pith(p)=p
	enddo
	poff=.true.
	
	call rcurve(mxdat,mxrec,nc(1),x1,xrc,idate,pith,poff,flag)
	mxage=maximum(age,nc(1))
	
       call detrend1(xrc,f,wk2,wk3,mxage,mxage,eq,jdt,mxdat,mxrec,ndt)			!Hakan

	call plt1(mxdat,mxrec,xrc,f,1,mxage,flag)	
		
	write(*,'(/1x,''Try another ?'',/1x,''Enter [1] for yes [2] for no :> '',$)')
	read(*,*)p
	if(p.eq.1)goto 40
	
	
c	
	write(*,*)
	write(*,'(13x,''Normal termination of program'')')
	write(*,'(13x,i4,'' files '')')fnum
	write(*,'(13x,i10,'' records '')')nc(1)
	write(*,*)
	stop	
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine plt1(mxsl,rmax,mey,ycal,fycaly,lycaly,flag)
c
c	argument
	integer mxsl,rmax
	integer fycaly,lycaly
	real*4 mey(mxsl),ycal(mxsl)
c	local
	integer j,k,i
	real*4 miny, maxy
	real*4 wrk1(500)
	character*20 tit
	character ans*1,filename*9,dev*80
	logical ok
c
       dev='Mac1';ok=.true.
c
10	if(ok)call plsdev(dev)
c
c	
	miny=9999.;maxy=-9999.
	do k=fycaly,lycaly
	    if(mey(k).ne.flag.and.mey(k).gt.maxy)maxy=mey(k)
	    if(mey(k).ne.flag.and.mey(k).lt.miny)miny=mey(k)
	    wrk1(k)=real(k)
	enddo
c
c initialize plplot
	call plinit()
c set lable color
	call plcol(15)   !color box and lables
c define the viewport one column by one rows
	call plssub(1,1)
	call plenv(real(fycaly),real(lycaly),miny,maxy,0,2)
	call pllab('BioYear','Y','Mean Curve and fit')
c	step=-12.
	call plcol(1)   !color box and lables
c	write(tit(1:20),'(i5)')nboot
c	call plmtex('rv',step,1.07,0.0,tit)
c	step=-6.
c	call plmtex('rv',step,1.07,0.0,'Reps')
	call plline(lycaly,wrk1,mey)
	call plcol(2)   !color box and lables
	call plline(lycaly,wrk1,ycal)
	call plend()
c	
       ans='C '
	dev='Mac1'
       write(*,'(/1x,''|************* plot options menu ***************|'')')
       write(*,'(1x,''|-----------------------------------------------|'')')
       write(*,'(1x,''|       [s]ave   <[c]ontinue>   [r]eplot        |'')') 
       write(*,'(1x,''|_______________________________________________|'')')
       write(*,'(/1x,'' enter option number > '',$)')
       read(*,'(a1)')ans
       if((ans.eq.'S').or.(ans.eq.'s'))then
         dev='psc';ok=.false.;filename='fig11.ps'
	  call plsdev(dev)
         call plsfnam(filename)
         goto 10
       elseif((ans.eq.'R').or.(ans.eq.'r'))then
         ok=.true.
         goto 10
       elseif((ans.eq.'C').or.(ans.eq.'c').or.(ans.eq.' ').or.
     *                             (ans.eq.'Y').or.(ans.eq.'y'))then
         return
	endif
c
	return
	end
c

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine rcurve(mxdat,mxrec,n,x1,xrc,idate,pith,poff,flag)
c compute a simple age aligned population mean curve
c argument decs
	integer mxdat,mxrec,n
	real*4 flag
	integer pith(mxdat),idate(mxrec,3)
	real*4 x1(mxdat,mxrec),xrc(mxdat)
	logical poff
c locals decs
	integer mxd, p,j,k, mxlyr
	parameter (mxd=500)    !must equal mxdat
	real*4  num(mxd),sum(mxd)  
c	
	mxlyr=-9999
	do j=1,mxd
	  sum(j)=0.0;num(j)=0.0
	enddo
	do j=1,n
	  p=pith(j)
	  do k=1,idate(j,3)
	     sum(p)=sum(p)+x1(k,j)
	     num(p)=num(p)+1.
	     p=p+1
	  enddo
	  if(idate(j,3).gt.mxlyr)mxlyr=idate(j,3)
	enddo
c	
	do p=1,mxlyr
	  if(num(p).ne.0)then
	    xrc(p)=sum(p)/num(p)
	  else
	    xrc(p)=flag
	  endif
c	  write(*,*)p,xrc(p)
	enddo
	return
	end
	
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr detrend1
*
      subroutine detrend1(x,f,wk2,wk3,n1,n,eq,jdt,mxy,mxs,ndt)			!Hakan
      parameter (i1=1,i2=2,i3=3,i4=4,i5=5,i6=6,i10=10,i21=21,iu11=11)
	  parameter (nemax=2500,ndmax=i5,ndmaxp1=i6,r0=0e0,r1=1e0,nmax=10000)
	  real*4 x(mxy),f(mxy),wk2(mxy),wk3(mxy),eq(7)
      real*4 pam(ndmaxp1,ndmax),yam(ndmaxp1),v(ndmax)
      real*4 span,alpha,w(nmax),sc(nmax,7)
      do i=1,7
        eq(i)=0.0
      enddo
      ls=0
      ndif=0
      if(jdt.eq.-5)then
        do i=1,n
          wk2(i)=float(i)
          wk3(i)=x(i)
          w(i)=1.0
        enddo
        span=0.0
        if(n.lt.40)span=0.3
        iper=1
        alpha=float(ndt)
        call supsmu(n,wk2,wk3,w,iper,span,alpha,f,sc)
      elseif(jdt.eq.-4)then
        do i=1,n
          wk2(i)=x(i)
        enddo
        rdt=float(ndt)/100.
c        rdt=float(ndt)
        call lowess(wk2,f,n,rdt)
      elseif(jdt.eq.-3)then
        do i=1,n
          wk2(i)=x(i)
        enddo
c        do i=2,n
c          wk2(i-1)=x(i)
c        enddo
        call msmoor(wk2,f,wk3,n,ndt,iter,ier,mxy)
      elseif(jdt.eq.-1)then
        do i=2,n
          wk2(i-1)=x(i)
        enddo
        call xmean(x(2),n-1,xm1)
        do i=2,n
          wk3(i-1)=x(i)-x(i-1)						                !first-difference
        enddo
        x(1)=-9.99
        n=n-1
        do i=1,n
          x(i+1)=wk3(i)
        enddo
c        call xmean(x,n,xm)
c        sum=0.0
c        do i=1,n
c          sum=sum+xm
c          f(i)=sum
c        enddo
c        call xmean(f,n,xm2)
c        do i=1,n
c          f(i)=f(i)-xm2+xm1
c        enddo
        ls=n/3
        call spline(n,wk2,f,p,ls)
        do i=1,n
          wk2(i)=wk2(i)/f(i)
        enddo
        call moment(x(2),n,ave1,adev,sdev1,var,skew,curt,sens,ac)
        call moment(wk2,n,ave2,adev,sdev2,var,skew,curt,sens,ac)
        do i=1,n
          x(i+1)=(x(i+1)-ave1)/sdev1
          x(i+1)=x(i+1)*sdev2+ave2
        enddo
      elseif(jdt.eq.1.or.jdt.eq.2)then
        do i=1,7
          eq(i)=0.0
        enddo
c        call curve(n,x,f,eq)				!negative exponential curve
        call curve(n1,n,x,f,eq)				!Hakan
        if(eq(2).le.0.or.eq(5).lt.0)then
          do i=1,7
            eq(i)=0.0
          enddo
c          call trend(n,x,f,eq(3),eq(5))
          call trend(n1,n,x,f,eq(3),eq(5))		!Hakan
          if(jdt.eq.1)then
            jdt=4
          elseif(jdt.eq.2.and.eq(3).le.0)then
            jdt=5
          else
            do i=1,7
              eq(i)=0.0
            enddo
            jdt=5
            call xmean(x,n,eq(5))
c            do i=1,n
            do i=1,n1
              f(i)=eq(5)					!Hakan
            enddo
          endif
        endif
      elseif(jdt.eq.3)then
        iopt=1
        do i=1,n
          wk2(i)=i
        enddo
	 nd=iopt+i2
        pam(1,1)=0.0
        pam(1,2)=0.0
        pam(1,3)=0.0
        pam(1,4)=0.0
        pam(1,5)=0.0
c        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n,wk2,x,f,sumsq)
        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n1,n,wk2,x,f,sumsq)		!Hakan
        slop=pam(i1,i2)
        eq(1)=pam(i1,i1)
        eq(2)=abs(pam(i1,i2))
        eq(3)=r0
        eq(5)=pam(i1,i3)
        eq(6)=r0
        if(slop.ge.0)then
          jdt=4
          do i=1,7
            eq(i)=0.0
          enddo
c          call trend(n,x,f,eq(3),eq(5))
          call trend(n1,n,x,f,eq(3),eq(5))				!Hakan
        endif
      elseif(jdt.eq.4.or.jdt.eq.5)then
c        call trend(n,x,f,eq(3),eq(5))				        !linear regression
        call trend(n1,n,x,f,eq(3),eq(5))				 	!Hakan
        if(jdt.eq.5.and.eq(3).gt.0)then
          jdt=6
          eq(3)=0.0
          eq(5)=0.0
          call xmean(x,n,eq(5))						            !arithmetic mean
c          do i=1,n
          do i=1,n1								!Hakan
            f(i)=eq(5)
          enddo 
        endif       
      elseif(jdt.eq.6)then
        call xmean(x,n,eq(5))						              !arithmetic mean
c        do i=1,n
        do i=1,n1								!Hakan
          f(i)=eq(5)
        enddo
      elseif(jdt.eq.7)then
        do i=1,n
          wk2(i)=i
        enddo
c        call hughdi(wk2,x,n,f,eq,flag)           !hugershoff growth curve
        call hughdi(wk2,x,n1,n,f,eq,flag)         !Hakan
        do i=1,n
          wk2(i)=i
          wk3(i)=f(i)
        enddo
        iopt=3
	       nd=iopt+i2
        pam(1,1)=eq(1)
        pam(1,2)=eq(2)
        pam(1,3)=eq(3)
        pam(1,4)=0.0
        pam(1,5)=0.0
c        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n,wk2,x,f,sumsq)
        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n1,n,wk2,x,f,sumsq)		!Hakan
        eq(1)=pam(i1,i1)
        eq(2)=pam(i1,i2)
        eq(3)=abs(pam(i1,i3))
        eq(5)=pam(i1,i5)
        eq(6)=pam(i1,i4)
c        if(eq(2).le.0)then
c          jdt=3
c          do i=1,7
c            eq(i)=0.0
c          enddo
c          call trend(n,x,f,eq(3),eq(5))
c        endif
        eq(7)=jdt
      elseif(jdt.eq.8)then
c        call gexp(x,n,f,wk2,wk3,eq,flag)         !general exponential curve
        call gexp(x,n1,n,f,wk2,wk3,eq,flag)       !Hakan
        do i=1,n
          wk2(i)=i
          wk3(i)=f(i)
        enddo
        iopt=2
	       nd=iopt+i2
        pam(1,1)=eq(1)
        pam(1,2)=eq(2)
        pam(1,3)=0.0
        pam(1,4)=0.0
        pam(1,5)=0.0
c        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n,wk2,x,f,sumsq)
        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n1,n,wk2,x,f,sumsq)		!Hakan
        eq(1)=pam(i1,i1)
        eq(2)=r1
        eq(3)=abs(pam(i1,i2))
        eq(5)=pam(i1,i4)
        eq(6)=pam(i1,i3)
        if(eq(2).le.0)then
          jdt=3
          do i=1,7
            eq(i)=0.0
          enddo
c          call trend(n,x,f,eq(3),eq(5))
          call trend(n1,n,x,f,eq(3),eq(5))			!Hakan
        endif
        eq(7)=jdt
      else
        if(jdt.gt.9)then
          ls=jdt											                      !fixed 50% frequency response
        elseif(jdt.lt.-9)then
          ls=ifix(abs((float(jdt)/100.0)*float(n)))			!%n frequency response
        endif
        do i=1,n
          wk2(i)=1.0
        enddo
ccc        call spl2w(x,n,f,wk2,p,ls,mxy)
ccc        call spl2w(mxy,n,x,wrk2,pp,f,ls)
        call spline(n,x,f,p,ls)						            !cubic spline filter
      endif
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      subroutine gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,ne,x,y,yfit,sumsq)
      subroutine gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n1,ne,x,y,yfit,sumsq)	!Hakan
c uses the nr subroutine amoeba to find useful values for the time and
c ringwidth offsets t and k in y = a*exp(b*ln(x+t)-c*(x+t)) + k as well
c as a, b, and c.
c
c nd is the number of amoeba dimensions. it must be ≤ ndmax.
c pam contains the nd+1 vectors of vertex coordinants of the amoeba.
c   they are initialized below to order of magnitude and usual sign.
c yam contains the rss for each vertex.
c iter is the number of steps the amoeba needed to find the solution.
c   itmax in amoeba is currently set to ????.
c if iopt = 1, an offset gamma with b=0 will be fit.
c  "   "  " 2,  "    "     "    "   b=1   "  "   " .
c           3,  "    "     "              "  "   " .
c ne is the number of (x,y) data pairs. it must be ≤ nemax.
c y is the (ringwidth) data to be fit.
c x is a vector of the same length as y. it will be loaded with the
c   numbers 1 thru ne below. not the years!
c yfit is the fitted function,
c sumsq it's residual sum of squares (rss).
c
c
      parameter (i1=1,i2=2,i3=3,i4=4,i5=5,i6=6)
      real*4    pam(ndmaxp1,ndmax),yam(ndmaxp1),v(ndmax)
	     real*4    pami(i3,i6,i5)
	     parameter (r0=0e0,rh=-1e-2,rt=1e-1,r1=1e0)
      parameter (ftol=1e-5)
	     real*4    x(nemax),y(nemax),yfit(nemax)
c234567890         0         0         0         0         0         012
	     data pami/r0,r0,r0,r1,r1,r1,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,
     *          r0,r0,r0,r0,r0,r0,rh,rh,r1,r0,r0,r0,r0,r0,r0,r0,r0,r0,
     *          r0,r0,r0,r0,r0,r0,r0,r0,r0,rt,r1,rh,r0,r0,r0,r0,r0,r0,
     *          r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,rt,r1,r0,r0,r0,
     *          r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,rt/
c
c initialize pam and yam.
c the first vertex has been passed in pam as input.
c replace the original(0) vertex with this linearized solution.
      do ic=i1,i5
	       v(ic)=pam(i1,ic)
	     enddo
	     yam(i1)=gamoff(v,ndmax,nemax,ne,x,y,iopt)
      do ir=i2,nd+i1
	       do ic=i1,nd
		        pamu=pami(iopt,ir,ic)
		        pam(ir,ic)=pamu
		        v(ic)=pamu
		      enddo
	       yam(ir)=gamoff(v,ndmax,nemax,ne,x,y,iopt)
	     enddo
c
c call amoeba.
    2 call amoeba(pam,yam,ndmaxp1,ndmax,nd,iter,ftol,
     *            iopt,nemax,ne,x,y)
	     sumsq=yam(i1)
c
c get the fitted function.
	     go to (21,22,23) iopt
   21 a=pam(i1,i1)
	     b=r0
      c=pam(i1,i2)
      t=r0
      rk=pam(i1,i3)
	     go to 9
   22 a=pam(i1,i1)
	     b=r1
	     c=pam(i1,i2)
	     t=pam(i1,i3)
	     rk=pam(i1,i4)
	     go to 9
   23 a=pam(i1,i1)
   	  b=pam(i1,i2)
   	  c=pam(i1,i3)
   	  t=pam(i1,i4)
   	  rk=pam(i1,i5)
c    9 do ie=i1,ne
    9 do ie=i1,n1						!Hakan
	       time=float(ie)+t
	       yfit(ie)=a*exp(b*alog(time)+c*time)+rk
	     enddo
	     return
	     end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	function gamoff(v,ndmax,nemax,ne,x,y,iopt)
c will compute y = a*exp(b*ln(x+t)-c*(x+t)) + k for x=1,ne.
c
      parameter (i1=1,i2=2,i3=3,i4=4,i5=5)
	     real*4    v(ndmax),x(nemax),y(nemax)
	     real*8    ad,bd,cd,td,rkd,sumsq,xd,tmd,d0,d1,yd,difd
	     parameter (d0=0d0,d1=1d0)
c
c dimensions...
      if(ne.gt.nemax) then
	       write(*,'(/1x,''in gamoff ne > nemax: '',2i10)') ne,nemax
		      pause'bye.'
		      stop
	     endif
c
c compute the residuals.
	     if(iopt.eq.i1) then
	       ad=dble(v(i1))
	       bd=d0
		      cd=dble(v(i2))
		      td=d0
		      rkd=dble(v(i3))
	       go to 8
	     endif
	     if(iopt.eq.i2) then
	       ad=dble(v(i1))
	       bd=d1
		      cd=dble(v(i2))
		      td=dble(v(i3))
		      rkd=dble(v(i4))
	       go to 8
	     endif
	     if(iopt.eq.i3) then
	       ad=dble(v(i1))
	       bd=dble(v(i2))
		      cd=dble(v(i3))
		      td=dble(v(i4))
		      rkd=dble(v(i5))
	       go to 8
	     endif
    8 sumsq=d0
	     do ie=i1,ne
	       xd=dble(x(ie))
		      tmd=xd+td
		      if(tmd.lt.d0) then
		        write(*,'(/1x,''in gamoff tmd ≤ 0 at xd.'',2d17.10)') tmd,xd
c23456789012345678901234567890123456789012345678901234567890123456789012
          pause'this is a fatal condition.'
		        stop
		      endif
		      yd=ad*dexp(bd*dlog(tmd)+cd*tmd)+rkd
		      difd=yd-dble(y(ie))
		      sumsq=sumsq+difd*difd
	     enddo
	     gamoff=sngl(sumsq)
	     return
	     end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine amoeba(p,y,mp,np,ndim,iter,ftol,iopt,
     *                 nemax,ne,x,ye)
      integer iter,mp,ndim,np,nmax,itmax
      real ftol,p(mp,np),y(mp)
	     real*4    x(nemax),ye(nemax)
      parameter (nmax=5,itmax=2000)
cu    uses ametry,gamoff
      integer i,ihi,ilo,inhi,j,m,n
      real rtol,sum,swap,ysave,ytry,psum(nmax)
      iter=0
1     do 12 n=1,ndim
        sum=0.
        do 11 m=1,ndim+1
          sum=sum+p(m,n)
11      continue
        psum(n)=sum
12    continue
2     ilo=1
      if (y(1).gt.y(2)) then
        ihi=1
        inhi=2
      else
        ihi=2
        inhi=1
      endif
      do 13 i=1,ndim+1
        if(y(i).le.y(ilo)) ilo=i
        if(y(i).gt.y(ihi)) then
          inhi=ihi
          ihi=i
        else if(y(i).gt.y(inhi)) then
          if(i.ne.ihi) inhi=i
        endif
13    continue
      rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo)))
      if (rtol.lt.ftol) then
        swap=y(1)
        y(1)=y(ilo)
        y(ilo)=swap
        do 14 n=1,ndim
          swap=p(1,n)
          p(1,n)=p(ilo,n)
          p(ilo,n)=swap
14      continue
        return
      endif
      if (iter.ge.itmax) then
        pause 'itmax exceeded in amoeba'
        return
      endif
      iter=iter+2
c	  write(*,'(1x,''iter,y: '',i5,3e14.7)') iter,(y(iy),iy=1,3)
c	  write(*,'(1x,''p: '')')
c	  do iv=1,3
c	    write(*,'(1x,2e14.7)') (p(iv,itk),itk=1,2)
c	  enddo
c	  write(*,'(1x,''a,b,c: '',3e14.7)') cc(1),cc(2),cc(3)
c234567890         0         0         0         0         0         012
      ytry=ametry(p,y,psum,mp,np,ndim,ihi,-1.0,iopt,
     *            nemax,ne,x,ye)
      if (ytry.le.y(ilo)) then
        ytry=ametry(p,y,psum,mp,np,ndim,ihi,2.0,iopt,
     *              nemax,ne,x,ye)
      else if (ytry.ge.y(inhi)) then
        ysave=y(ihi)
        ytry=ametry(p,y,psum,mp,np,ndim,ihi,0.5,iopt,
     *              nemax,ne,x,ye)
        if (ytry.ge.ysave) then
          do 16 i=1,ndim+1
            if(i.ne.ilo)then
              do 15 j=1,ndim
                psum(j)=0.5*(p(i,j)+p(ilo,j))
                p(i,j)=psum(j)
15            continue
			  itag=8
c234567890         0         0         0         0         0         012
              y(i)=gamoff(psum,np,nemax,ne,x,ye,iopt)
            endif
16        continue
          iter=iter+ndim
          goto 1
        endif
      else
        iter=iter-1
      endif
      goto 2
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      function ametry(p,y,psum,mp,np,ndim,ihi,fac,iopt,
     *                nemax,ne,xe,ye)
      integer ihi,mp,ndim,np,nmax
      real fac,p(mp,np),psum(np),y(mp)
	     real*4    xe(nemax),ye(nemax)
      parameter (nmax=5)
cu    uses gamoff
      integer j
      real fac1,fac2,ytry,ptry(nmax)
      fac1=(1.-fac)/ndim
      fac2=fac1-fac
      do 11 j=1,ndim
        ptry(j)=psum(j)*fac1-p(ihi,j)*fac2
11    continue
      i1=1
	     ndimm1=ndim-i1
	     r0=0e0
      if(iopt.gt.i1) then
		      if(ptry(ndimm1).le.r0) ptry(ndimm1)=r0
	     endif
	     itag=9
      ytry = gamoff(ptry,np,nemax,ne,xe,ye,iopt)
      if (ytry.lt.y(ihi)) then
        y(ihi)=ytry
        do 12 j=1,ndim
          psum(j)=psum(j)-p(ihi,j)+ptry(j)
          p(ihi,j)=ptry(j)
12      continue
      endif
      ametry=ytry
      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr curve: compute exponential curve coefficients
*
*  programmed by fritts, huzar, bottorf and ray
*  revised oct 1968 by v bockman
*  revised     1983 by e r cook
*  revised by r l holmes    feb 1984
*  revised by r l holmes  5 feb 1985
*
*     tree-ring laboratory
*     university of arizona
*     tucson, arizona  85721  usa
*
***********************************************************
*
* parameters passed to routine:
*
* n:  number of data in z
* z:  array of data to which to fit curve
* f:  working array of length n
* eq: array of exponential curve coefficients computed here
*
* exponential curve equation computed:
*
*     y = ah * exp(-bh * t) + eh
*        eq(1)     eq(2)     eq(5)
* eq(6): iterations required to converge
* eq(7): sum of array z
*
**************************************************
*
c      subroutine curve(n,z,f,eq)
      subroutine curve(n1,n,z,f,eq)		!Hakan
      real*4 z(n),f(n),eq(7)
      real*8 a,b,d,e,g,eh,ah,bh(3),rss(3),yf,yl,rn,dlt,ab,arg,zz,dinc
      data dinc/.0000005d0/
      nit=0
      if(n .lt. 20)then
      write(*,'(''>> sbr curve:  only'',i3,
     +'' values, no exponential curve fit'')')n
      goto 99
      endif
      rn=dble(n)
      yf=0.d0
      do 5 i=1,10
    5 yf=yf+dble(z(i))
      yf=yf/10.d0
      yl=0.d0
      do 8 i=n-9,n
    8 yl=yl+dble(z(i))
      yl=yl/10.d0
      bh(2)=(dlog(yf)-dlog(yl))/(rn-10.d0)
      if(bh(2) .le. 0.)bh(2)=1.d-8
      dlt=.0005d0
      mu=0
      mo=0
*
* iteration loop
      do 13 it=1,34
      nit=it
      bh(1)=bh(2)-dlt
      bh(3)=bh(2)+dlt
      do 28 k=1,3
      if(k .eq. 2 .and. mu .gt. 0)goto 28
      ab=abs(bh(k))
      do 19 i=1,n
      arg=ab*dfloat(i)
      if(arg .gt. 112.8d0)goto 99
   19 f(i)=dexp(-bh(k)*dfloat(i))
      a=dble(n)
      b=0.d0
      d=0.d0
      e=0.d0
      g=0.d0
      do 21 i=1,n
      b=b+dble(f(i))
      d=d+dble(z(i))
      e=e+dble(f(i))*dble(f(i))
   21 g=g+dble(z(i))*dble(f(i))
      if((b .gt. 1.d15 .or. e .gt. 1.d15) .or.
     +(b .lt. 1.d-15 .or. e .lt. 1.d-15))goto 99
      eh=(d*e-b*g)/(a*e-b*b)
      ah=(d-a*eh)/b
      rss(k)=0.d0
      do 25 i=1,n
      zz=dble(z(i))-eh-ah*dble(f(i))
   25 rss(k)=rss(k)+zz*zz
   28 continue
      mu=1
      lo=0
      do 30 k=1,3,2
      if(rss(k) .lt. rss(2))then
      lo=-1
      bh(2)=bh(k)
      rss(2)=rss(k)
      endif
   30 continue
      if(lo .ge. 0)mo=1
      if(mo .le. 0)then
      dlt=dlt*2.d0
      else
      dlt=dlt*.5d0
      if(dlt .le. dble(dinc))goto 15
      endif
   13 continue
      write(*,'(''>> sbr curve: over 34 iterations'')')
      goto 99
   15 if(eh.gt.9999.d0 .or. ah.le.0.d0 .or. bh(2).lt.0.d0)bh(2)=0.d0
      eq(1)=sngl(ah)
      eq(2)=bh(2)
      eq(5)=sngl(eh)
      eq(6)=nit
      eq(7)=sngl(d)
*
* calculate curve values
c      do i=1,n
      do i=1,n1							!Hakan
        f(i)=sngl(ah*dexp(bh(2)*dfloat(i)*(-1.d0))+eh)
      enddo
      return
*
* flag values if coefficients not computed
   99 f(1)=-999.
      eq(2)=-999.
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr trend: fits trend line to series
*
*  programmed oct 1982 by
*
*     richard  l  holmes
*     tree-ring laboratory
*     university of arizona
*     tucson, arizona  85721  usa
*
*  modified by r l holmes 15 mar 1985
*  modified by e r cook 25 aug 1985 to double precision
*
********************************************************************
*
*  n:  number of values in series
*  y:  time series array
*  f:  trend line values calculated by this routine
*  sl: slope of trend line
*  yi: y-axis intercept
*
********************************************************************
*
c      subroutine trend(n,y,f,sl,yi)
      subroutine trend(n1,n,y,f,sl,yi)			!Hakan
      implicit real*8 (a-h,o-z)
      real*4 y(n),f(n),sl,yi
      an=dfloat(n)
      t=0.d0
      w=0.d0
      x=0.d0
      z=0.d0
      do i=1,n
        ai=dfloat(i)
        z=z+ai
        t=t+ai*ai
        x=x+dble(y(i))*ai
        w=w+dble(y(i))
      enddo
      ymean=w/an
      dsl=(x-w*z/an)/(t-z*z/an)
      dyi=ymean-dsl*z/an
c      do i=1,n
      do i=1,n1						!Hakan
        ai=dfloat(i)
        f(i)=sngl(dsl*ai+dyi)
      enddo
      sl=sngl(dsl)
      yi=sngl(dyi)
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr hughdi
*
c      subroutine hughdi(x,y,ntot,xy,coef,flag)
      subroutine hughdi(x,y,n1,ntot,xy,coef,flag)		!Hakan
c
c     x     --  predictor variable -- time t=1,ntot
c     y     --  tree-ring series
c     ntot  --  number of years
c     xy    --  hugershoff growth curve estimates
c     coef  --  vector of hugershoff coefficients
c     flag  --  missing value flag
c
c     ausgleich mit hugershoff-funktion
c     schaetzen der parameter mit direkter berechnung (h.p. holzhauser)
c     a, cm, ck    im modell 
c     y(i) = a * t(i)**cm * exp(ck*t(i))   fuer i=1....ntot
c     fuer die berechnung muessen folgender summenregister vom 
c     logarithmierten x und y und unlogarithmierten x (t) vorbereitet werden.
c     sumt, sumt2, sumx, sumx2, sumy, sumyt, sumxy und sumxt.
c     regression ueber hilfsfunktion 
c     z = ln(y) = ln(a) + m*ln(t) - k*t 
c     umgesetzt in z = a' + m*t' - k*t
c     nrow entsprechend dimensionierung hauptprogramm
c     ncoef muss 3 sein
c
c     september 1989 / nos/ve /  braeker/nogler
c     pc-version 2.1 / juni 1992 / p. nogler
c     modified by e.r. cook oct 23, 1996
c
      implicit real*8 (a-h,o-z)
      real*4 x(ntot),y(ntot),xy(ntot),coef(7),flag
      real*8 a,ck,cm,se,sf,cc,dcoef(3)
      integer ntot,ncoef,nred
      logical lber
      lber=.true.
      cc = -999.0d0
c
c     beginn hughdi
c
      if(.not.lber) goto 500 
      nred  = 0 
      sumt  = 0.0d0
      sumx  = 0.0d0
      sumx2 = 0.0d0
      sumt2 = 0.0d0
      sumxt = 0.0d0
      sumy  = 0.0d0
      sumxy = 0.0d0
      sumyt = 0.0d0
c
      do 200 k = 1,ntot
        xk=dble(x(k))
        yk=dble(y(k))
        if(yk.eq.flag) goto 200
        if(yk.le.0.0d0) goto 200
        if(xk.le.0.0d0) goto 200
        nred = nred+1 
        xl = dlog(xk)
        yl = dlog(yk)
        sumt  = sumt  + xk
        sumt2 = sumt2 + xk*xk
        sumx  = sumx  + xl
        sumx2 = sumx2 + xl*xl
        sumxt = sumxt + xl*xk
        sumy  = sumy  + yl
        sumxy = sumxy + xl*yl
        sumyt = sumyt + yl*xk
  200 continue
c 
      rn1 = (sumx2 * sumt2 - sumxt**2) 
      rn2 = (sumt * sumxt - sumx * sumt2) 
      rn3 = (sumx * sumxt - sumt * sumx2)
      rn  = nred * rn1 + sumx * rn2 + sumt * rn3
      rza = sumy * rn1 + sumxy * rn2 + sumyt * rn3
      a   = rza / rn
      rz2 = (sumx * sumt - nred * sumxt)
      rzm = sumy * rn2 + sumxy * (nred * sumt2 - sumt**2) + sumyt * rz2
      cm  = rzm / rn
      rz1 = (sumt * sumx2 - sumx * sumxt)
      rz2 = (nred * sumxt - sumx * sumt)
      rzk = sumy * rz1 + sumxy * rz2 + sumyt * (sumx**2 - nred * sumx2)
      ck  = rzk / rn
c
      ncoef = 3
      dcoef(1) = dexp(a)
      dcoef(2) = cm
      dcoef(3) = -ck
      coef(1) = sngl(dexp(a))
      coef(2) = sngl(cm)  
      coef(3) = sngl(-ck) 
ccc      if(coef(2).le.0.0)then
ccc        coef(2)=-999.0
ccc        xy(1)=-999.0
ccc        goto 700
ccc      endif
c 
c     coef-korrekturfaktor berechnen ueber abweichungen
c
      se = 0.0d0
      sf = 0.0d0
c      do k = 1,ntot 
      do k = 1,n1							!Hakan 
        xk = dble(x(k))
        yk = dble(y(k))
        dxy = dcoef(1) * xk**dcoef(2) * dexp(dcoef(3)*xk)
        xy(k) = sngl(dxy)
        if(yk.ne.flag)then 
          se = se+(yk-dxy) 
          sf = sf+dxy 
        endif
      enddo
      cc = se/sf+1.0d0
      if(cc.eq.0.0d0) goto 700 
      dcoef(1) = dcoef(1)*cc
      coef(1) = sngl(dcoef(1))
c      do k = 1,ntot 
      do k = 1,n1							!Hakan 
        xy(k) = xy(k)*cc 
      enddo
      goto 700
  500 continue
c     berechnen der ausgleichskurve und xy in kolonne 1 setzen
c      do k = 1,ntot 
      do k = 1,n1							!Hakan
        xk = dble(x(k))
        dxy = dcoef(1) * xk**dcoef(2) * dexp(dcoef(3)*xk)
        xy(k) = sngl(dxy) 
      enddo
  700 return
      end 
*sbr gexp
*
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine gexp(y,nt,n,f,wk1,wk2,coef,flag)		!Hakan
c
c     y     --  tree-ring series
c     n     --  number of years
c     f     --  general exponential curve estimates
c     wk1   --  work vector
c     wk2   --  work vector
c     coef  --  vector of hugershoff coefficients
c     flag  --  missing value flag
c
      real*4 y(n),f(n),wk1(n),wk2(n),coef(7),flag
      real*8 a,b,r,ax,ay,dn1,xt,yt,sxx,syy,sxy,sgy,sgx,yxr,ri
      if(n.le.0)then
        write(*,'(/1x,''*** n .le. 0 in sbr gexp ***''/)')
        pause
        stop
      endif
      a=0.0d0
      b=0.0d0
      n1=0
c
c *** linearize the problem
c
      do i=1,n
        if(y(i).gt.0.and.y(i).ne.flag)then
          n1=n1+1
          wk1(n1)=float(i)
          wk2(n1)=alog(y(i))-alog(wk1(i))
        endif
      enddo
      if(n1.le.0)then
        write(*,'(/1x,''*** n1 .le. 0 in sbr gexp ***''/)')
        pause
        stop
      endif
c
c *** now compute simple linear regression
c
      dn1=dfloat(n1)
      ax=0.0d0
      ay=0.0d0
      do j=1,n1
        ax=ax+dble(wk1(j))
        ay=ay+dble(wk2(j))
      enddo
      ax=ax/dn1
      ay=ay/dn1
      sxx=0.0d0
      syy=0.0d0
      sxy=0.0d0
      do j=1,n1
        xt=dble(wk1(j))-ax
        yt=dble(wk2(j))-ay
        sxx=sxx+xt**2
        syy=syy+yt**2
        sxy=sxy+xt*yt
      enddo
      sgy=dsqrt(syy/(dn1-1.0d0))
      sgx=dsqrt(sxx/(dn1-1.0d0))
      yxr=sgy/sgx
      r=sxy/dsqrt(sxx*syy)
      b=yxr*r
      a=ay-b*ax
c
c *** unlinearize the fitted model and calculate growth curve values
c
      a=dexp(a)
c      do i=1,n
      do i=1,nt					!Hakan
        ri=dble(i)
        f(i)=sngl(a*ri*dexp(b*ri))
      enddo
      coef(1)=sngl(a)
      coef(2)=-sngl(b)
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr moment
*
      subroutine moment(data,n,ave,adev,sdev,var,skew,curt,sens,ac)
      real*4 data(n)
      if(n.eq.0)then
        ave=-99.0
        adev=-99.0
        sdev=-99.0
        var=-99.0
        skew=-99.0
        curt=-99.0
        sens=-99.0
        ac=-99.0
        return
      endif
      rn=n
      s=0.0
      do j=1,n
        s=s+data(j)
      enddo
      ave=s/rn
      if(n.eq.1)then
        adev=-99.0
        sdev=-99.0
        var=-99.0
        skew=-99.0
        curt=-99.0
        sens=-99.0
        ac=-99.0
        return
      endif
      adev=0.0
      var=0.0
      skew=0.0
      curt=0.0
      sens=0.0
      do j=1,n
        s=data(j)-ave
        adev=adev+abs(s)
        p=s*s
        var=var+p
        p=p*s
        skew=skew+p
        p=p*s
        curt=curt+p
        if(j.lt.n)then
          rnum=data(j+1)-data(j)
          dnum=(data(j+1)+data(j))/2.0
          if(dnum.gt.0)then
            sens=sens+abs(rnum/dnum)
          endif
        endif
      enddo
      sens=sens/(rn-1.0)
      adev=adev/rn
      var=var/(rn-1.0)
      bvar=var*(rn-1.0)/rn
      if(var.gt.0.0)then
        sdev=sqrt(var)
        bsdev=sqrt(bvar)
        skew=skew/(rn*sdev**3)
        curt=curt/(rn*var**2)
        if(n.lt.3)then
          bfs=1.0
        else
          bfs=(sqrt(abs((rn*(rn-1.0))))/(rn-2.0))*(1.0+8.5/rn)
        endif
        if(n.lt.4)then
          bfk=1.0
        else
          bfk=(rn*rn*(rn+1.0))/((rn-1.0)*(rn-2.0)*(rn-3.0))
        endif
        skew=skew*bfs
        curt=curt*bfk
        do i=1,n
	         data(i)=(data(i)-ave)/bsdev
        enddo
        cxx=0.0
        do i=2,n
          cxx=cxx+data(i)*data(i-1)
        enddo
        ac=cxx/rn
        do i=1,n
          data(i)=data(i)*bsdev+ave
        enddo
      else
        sdev=-99.0
        skew=-99.0
        curt=-99.0
        ac=-99.0
      endif
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr spline: fits cubic spline curve to time series
*
*  derived from imsl library routines by edward r cook,
*  lamont-doherty geological observatory, and modified by
*  richard l holmes, tree-ring laboratory, university of arizona.
*
*  modified by r l holmes 10 jun 1988
*  modified by e r cook to full double precision 4 apr 1991
*
********************************************************************
*
* calculates a cubic spline function to model a time series.
* 'stiffness' of spline is specified by parameter 'ls'
*
* parameters:
*
* n:    number of values in time series
* y:    time series array to be modeled
* f:    array of cubic spline function values computed
* p:    value computed in this routine, used to calculate spline
* ls:   length (stiffness) of spline to be used to model series
* a:    working array of dimension (n,4)
*
********************************************************************
*
      subroutine spline(n,y,f,pp,ls)
      parameter(mxy=5000)
      implicit real*8 (a-h,o-z)
      dimension a(mxy,4),c1(4),c2(3)
      real*4 y(n),f(n),pp
      data c1 /1.d0,-4.d0,6.d0,-2.d0/
      data c2 /0.d0,.33333333333333d0,1.33333333333333d0/
      data pct,pi / .50d0, 3.1415926535897935d0 /
      if(n .lt. 4)then
      f(1)=-9998.
      return
      endif
      if(n.gt.mxy)then
        write(*,'(/t2,''*** series too long for spline subroutine ***''/
     *           t2,'' *** enlarge mxy parameter in subroutine ***''/)')
        stop
      endif
      nm2=n-2
      v=dfloat(ls)
      arg=(2.d0*pi)/v
      p=(6.d0*(dcos(arg)-1.d0)**2.d0)/(dcos(arg)+2.d0)
      pp=sngl(p)
ccc      p1=((1.d0/(1.d0-pct)-1.d0)*6.d0*(dcos(pi*2.d0/v)-1.d0)**2)/
ccc     +     (dcos(pi*2.d0/v)+2.d0)
ccc      p2=(6.d0*(dcos(arg)-1.d0)**2.d0)/((2.d0**0.5d0-1.d0)*
ccc     *		(dcos(arg)+2.d0))
      do 1 i=1,nm2
      do 1 j=1,3
      a(i,j)=c1(j)+p*c2(j)
    1 a(i,4)=dble(y(i))+c1(4)*dble(y(i+1))+dble(y(i+2))
      a(1,1)=c2(1)
      a(1,2)=c2(1)
      a(2,1)=c2(1)
      nc=2
* begin ludapb
      rn=dble(1.d0/(dfloat(nm2)*16.d0))
      d1=1.d0
      d2=0.d0
      ncp1=nc+1
      if(nc .ne. 0)then
* initialize zero elements
      do 5 i=1,nc
      do 5 j=i,nc
      k=ncp1-j
    5 a(i,k)=0.d0
* 'i' is row index of element being computed
* 'j' is column index of element being computed
* 'l' is row index of previously computed vector
*       being used to compute inner product
* 'm' is column index
      endif
      do 60 i=1,nm2
      imncp1=i-ncp1
      i1=max0(1,1-imncp1)
      do 60 j=i1,ncp1
      l=imncp1+j
      i2=ncp1-j
      sum=a(i,j)
      jm1=j-1
      if(jm1 .gt. 0)then
      do 25 k=1,jm1
      m=i2+k
   25 sum=sum-a(i,k)*a(l,m)
      endif
      if(j .eq. ncp1)then
      if(a(i,j)+sum*rn .le. a(i,j))then
      f(1)=-9999.
      write(*,'(''>> sbr spline: matrix not positive definite'')')
      return
      endif
      a(i,j)=1.d0/dsqrt(sum)
* update determinant
      d1=d1*sum
   35 if(dabs(d1) .gt. 1.d0)then
      d1=d1*.0625d0
      d2=d2+4.d0
      goto 35
      endif
   47 if(dabs(d1) .le. .0625d0)then
      d1=d1*16.d0
      d2=d2-4.d0
      goto 47
      else
      goto 60
      endif
      endif
      a(i,j)=sum*a(l,ncp1)
   60 continue
* end ludapb / begin luelpb
* solution ly = b
      nc1=nc+1
      iw=0
      l=0
      do 15 i=1,nm2
      sum=a(i,4)
      if(nc .gt. 0)then
      if(iw .ne. 0)then
      l=l+1
      if(l .gt. nc) l=nc
      k=nc1-l
      kl=i-l
      do 45 j=k,nc
      sum=sum-a(kl,4)*a(i,j)
   45 kl=kl+1
      else
      if(sum .ne. 0.d0)iw=1
      endif
      endif
   15 a(i,4)=sum*a(i,nc1)
* solution ux = y
      a(nm2,4)=a(nm2,4)*a(nm2,nc1)
      n1=nm2+1
      do 30 i=2,nm2
      k=n1-i
      sum=a(k,4)
      if(nc .gt. 0)then
      kl=k+1
      k1=min0(nm2,k+nc)
      l=1
      do 50 j=kl,k1
      sum=sum-a(j,4)*a(j,nc1-l)
   50 l=l+1
      endif
   30 a(k,4)=sum*a(k,nc1)
* end luelpb
      do 2 i=3,nm2
    2 f(i)=a(i-2,4)+c1(4)*a(i-1,4)+a(i,4)
      f(1)=a(1,4)
      f(2)=c1(4)*a(1,4)+a(2,4)
      f(n-1)=a(nm2-1,4)+c1(4)*a(nm2,4)
      f(n)=a(nm2,4)
      do 3 i=1,n
    3 f(i)=y(i)-f(i)
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr xmean
*
      subroutine xmean(x,n,xm)
      real*4 x(n)
      if(n.le.0)then
        write(*,'(/1x,''error in xmean: n=0'')')
        pause
        stop
      endif
      rn=n
      xm=0.0
      do i=1,n
        xm=xm+x(i)
      enddo
      xm=xm/rn
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr lowess
*
*  argument description:
* 
*       x = input; abscissas of the points on the scatterplot; 
*           the values in x must be ordered from smallest to largest.
*       y = input; ordinates of the points on the scatterplot.
*       n = input; dimension of x,y,ys,rw, and res.
*       f = input; specifies the amount of smoothing; f is
*           the fraction of points used to compute each
*           fitted value; as f increases the smoothed values
*           become smoother; choosing f in the range .2 to .8
*           usually results in a good fit; if you have no
*           idea which value to use, try f = .5.
*  nsteps = input; the number of iterations in the robust
*           fit; if nsteps = 0, the nonrobust fit is
*           returned; setting nsteps equal to 2 should serve
*           most purposes.
*   delta = input; nonnegative parameter which may be used
*           to save computations; if n is less than 100, set
*           delta equal to 0.0; if n is greater than 100 you
*           should find out how delta works by reading the
*           additional instructions section.
*      ys = output; fitted values; ys(i) is the fitted value
*           at x(i); to summarize the scatterplot, ys(i)
*           should be plotted against x(i).
*      rw = output; robustness weights; rw(i) is the weight
*           given to the point (x(i),y(i)); if nsteps = 0,
*           rw is not used.
*     res = output; residuals; res(i) = y(i)-ys(i).
* 
ccc      subroutine lowess(x,y,n,f,nsteps,delta,ys,rw,res)
      subroutine lowess(y,ys,n,f)
      parameter(mxy=5000)
      integer n
      integer nsteps
ccc      real x(n), y(n), f, delta, ys(n), rw(n)
ccc      real res(n)
      real x(mxy), y(n), f, delta, ys(n), rw(mxy)
      real res(mxy)
      integer nright, min0, max0, i, j, ifix
      integer iter, last, m1, m2, ns, nleft
      real abs, cut, cmad, r, d1, d2
      real c1, c9, alpha, denom, float
      logical ok
c
c *** ed's add in stuff
c
      nsteps=2
c	     if(n.le.100)then
        delta=0.0
c      else
c        range=x(n)-x(1)
c        delta=range/50.0
c      endif
      do i=1,n
        ri=i
        x(i)=ri
      enddo
c
c ***
c
      if (n .ge. 2) goto 1
         ys(1) = y(1)
         return
c at least two, at most n points
   1  ns = max(min(int(f*float(n)), n), 2)   !pjk
c   1  ns = max0(min0(ifix(f*float(n)), n), 2)
      iter = 1
         goto  3
   2     iter = iter+1
   3     if (iter .gt. nsteps+1) goto  22
c robustness iterations
         nleft = 1
         nright = ns
c index of prev estimated point
         last = 0
c index of current point
         i = 1
   4        if (nright .ge. n) goto  5
c move nleft, nright to right if radius decreases
               d1 = x(i)-x(nleft)
c if d1<=d2 with x(nright+1)==x(nright), lowest fixes
               d2 = x(nright+1)-x(i)
               if (d1 .le. d2) goto  5
c radius will not decrease by move right
               nleft = nleft+1
               nright = nright+1
               goto  4
c fitted value at x(i)
   5        call lowest(x, y, n, x(i), ys(i), nleft, nright, res, iter
     +     .gt. 1, rw, ok)
            if (.not. ok) ys(i) = y(i)
c all weights zero - copy over value (all rw==0)
            if (last .ge. i-1) goto 9
               denom = x(i)-x(last)
c skipped points -- interpolate
c non-zero - proof?
               j = last+1
                  goto  7
   6              j = j+1
   7              if (j .ge. i) goto  8
                  alpha = (x(j)-x(last))/denom
                  ys(j) = alpha*ys(i)+(1.0-alpha)*ys(last)
                  goto  6
   8           continue
c last point actually estimated
   9        last = i
c x coord of close points
            cut = x(last)+delta
            i = last+1
               goto  11
  10           i = i+1
  11           if (i .gt. n) goto  13
c find close points
               if (x(i) .gt. cut) goto  13
c i one beyond last pt within cut
               if (x(i) .ne. x(last)) goto 12
                  ys(i) = ys(last)
c exact match in x
                  last = i
  12           continue
               goto  10
c back 1 point so interpolation within delta, but always go forward
  13        i = max(last+1, i-1)
c  13        i = max0(last+1, i-1)
  14        if (last .lt. n) goto  4
c residuals
         do  15 i = 1, n
            res(i) = y(i)-ys(i)
  15        continue
         if (iter .gt. nsteps) goto  22
c compute robustness weights except last time
         do  16 i = 1, n
            rw(i) = abs(res(i))
  16        continue
         call sort(n,rw)
         m1 = n/2+1
         m2 = n-m1+1
c 6 median abs resid
         cmad = 3.0*(rw(m1)+rw(m2))
         c9 = .999*cmad
         c1 = .001*cmad
         do  21 i = 1, n
            r = abs(res(i))
            if (r .gt. c1) goto 17
               rw(i) = 1.
c near 0, avoid underflow
               goto  20
  17           if (r .le. c9) goto 18
                  rw(i) = 0.
c near 1, avoid underflow
                  goto  19
  18              rw(i) = (1.0-(r/cmad)**2)**2
  19        continue
  20        continue
  21        continue
         goto  2
  22  return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine lowest(x,y,n,xs,ys,nleft,nright,w,userw,rw,ok)
      parameter(mxy=5000)
      integer n
      integer nleft, nright
ccc      real x(n), y(n), xs, ys, w(n), rw(n)
      real x(mxy), y(n), xs, ys, w(n), rw(n)
      logical userw, ok
      integer nrt, j
      real abs, a, b, c, h, r
      real h1, sqrt, h9, range   !, amax1   -pjk
      range = x(n)-x(1)
      h = max(xs-x(nleft), x(nright)-xs)    !pjk
c      h = amax1(xs-x(nleft), x(nright)-xs)
      h9 = .999*h
      h1 = .001*h
c sum of weights
      a = 0.0
      j = nleft
         goto  2
   1     j = j+1
   2     if (j .gt. n) goto  7
c compute weights (pick up all ties on right)
         w(j) = 0.
         r = abs(x(j)-xs)
         if (r .gt. h9) goto 5
            if (r .le. h1) goto 3
               w(j) = (1.0-(r/h)**3)**3
c small enough for non-zero weight
               goto  4
   3           w(j) = 1.
   4        if (userw) w(j) = rw(j)*w(j)
            a = a+w(j)
            goto  6
   5        if (x(j) .gt. xs) goto  7
c get out at first zero wt on right
   6     continue
      
ccc         if(y(j).eq.-9.99)w(j)=0.0               !ed's bodge

         goto  1
c rightmost pt (may be greater than nright because of ties)
   7  nrt = j-1
      if (a .gt. 0.0) goto 8
         ok = .false.
         goto  16
   8     ok = .true.
c weighted least squares
         do  9 j = nleft, nrt
c make sum of w(j) == 1
            w(j) = w(j)/a
   9        continue
         if (h .le. 0.) goto 14
            a = 0.0
c use linear fit
            do  10 j = nleft, nrt
c weighted center of x values
               a = a+w(j)*x(j)
  10           continue
            b = xs-a
            c = 0.0
            do  11 j = nleft, nrt
               c = c+w(j)*(x(j)-a)**2
  11           continue
            if (sqrt(c) .le. .001*range) goto 13
               b = b/c
c points are spread out enough to compute slope
               do  12 j = nleft, nrt
                  w(j) = w(j)*(b*(x(j)-a)+1.0)
  12              continue
  13        continue
  14     ys = 0.0
         do  15 j = nleft, nrt
            ys = ys+w(j)*y(j)
  15        continue
  16  return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr msmoor
*
      subroutine msmoor(x,y,z,n,m,iter,ier,mxy)
c     repeatedly calls msmoo until either absum < small or iter > itmax.
c
      parameter (i0=0,i1=1,i2=2,i3=3)
      parameter (itmax=10,small=1e-6)                                                                               
      real*4  x(mxy),y(mxy),z(mxy)
c
c input error checks
      ier=i0
      if(n.lt.m.or.n.gt.mxy)ier=i1
      if(m.lt.i1.or.m.gt.mxy)ier=i2
      if(mod(m,i2).ne.i1)m=m-1
      if(mod(m,i2).ne.i1)ier=i3
      if(ier.ne.i0) then
        write(*,'(/1x,''input error '',i1,'' in msmoor.'')') ier
        pause
        return
      endif
c pass.
      iter=i0
   10 iter=iter+i1
      if(iter.gt.itmax) then
        ier=20
		      write(*,'(/1x,''*** maximum iterations reached in '',
     *    ''msmoor ***'')')
        return
      endif
      call msmoo(x,y,z,n,m,iter,ier,mxy)
      if(ier.ne.i0) then
        ier=ier+10
        return
      endif
      absum=0e0
      do i=i1,n                                                              
        absum = absum + abs(x(i) - y(i))
      enddo
c      write(*,'(/1x,i5,f10.7)') iter,absum
c      write(11,'(/1x,i5,f10.7)') iter,absum
      if(absum.lt.small) return
      do i=i1,n
        x(i)=y(i)
      enddo
      go to 10
      end 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr msmoo
*                                                                      
      subroutine msmoo(x,y,z,n,m,iter,ier,mxy)
c     will do median smoothing of a time series. n is the number of
c     data in the series, i.e. the series length. m is the filter
c     length which must be odd. x comes in,
c                               y goes out, and
c                               z is a little work vector.
c                               oh la di da!
c     at the ends medians of 1, 3, 5,..., m-2 data points are taken.
c
      parameter (i0=0,i1=1,i2=2)
      parameter (r0=0e0)
      real*4 x(mxy),y(mxy),z(mxy)
c
c
c input error checks
      ier=i0
      if(n.lt.m.or.n.gt.mxy) ier=i1
      if(m.lt.i1.or.m.gt.mxy) ier=i2
      if(mod(m,i2).ne.i1) ier=3
      if(ier.ne.i0) then
        write(*,'(/1x,''input error '',i1,'' in msmoo.'')') ier
        pause
        return
      endif
      if(m.eq.i1) then
        do i=i1,n
          y(i)=x(i)
        enddo
        return
      endif
      mm = (m-i1)/i2
c
c save the first point.
ccc      y(i1)=x(i1+1)
      y(i1)=x(i1)
c
c do from 2 thru mm.
      if(mm.gt.i1) then
        do i=i2,mm
          im=i2*i-i1
          do j=i1,im
            z(j) = x(j)
          enddo
ccc          call rank(z,im,mmax)
          call sort(im,z)
          y(i) = z(i)
        enddo
      endif
c
c do the middle.
      do i=mm+i1,n-mm
        jj=i-mm-i1
        do j=i-mm,i+mm
          z(j-jj) = x(j)
        enddo
ccc        call rank(z,m,mmax)
        call sort(m,z)
        y(i) = z(mm+i1)
      enddo
c
c do from n-mm+1 thru n-1.
      if(mm.gt.i1) then
        np1=n+i1
        do i=i2,mm
          im=i2*i-i1
          do j=i1,im
            z(j) = x(np1-j)
          enddo
ccc          call rank(z,im,mxy)
          call sort(im,z)
          y(np1-i)=z(i)
        enddo
      endif
c
c save the last point.
ccc      y(n)=x(n-1)
      y(n)=x(n)
c
      do i=i1,mm
        j=n-i+1
        y(i)=y(mm+1)
        y(j)=y(n-mm)
      enddo
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr trrw: read/write indices/measurements
*
*  programmed jun 1984 by
*
*     richard  l  holmes
*     tree-ring laboratory
*     university of arizona
*     tucson, arizona  85721  usa
*
*  modified by r l holmes 22 aug 1988
*
************************************************************
*
*  this routine reads or writes tree-ring
*  indices or measurements in standard formats used
*  in the international tree-ring data bank (itrdb).
*
*     call trrw(id,iyr,n,y,nc,ifl,iu,iuo,q)
*               ==       = ==            =        (= see below)
*
* arguments (9):
*
* id:    series identification
*           6 characters for indices, 8 for measurements; string may be
*           longer.  when q = 'ir', if the first six characters of
*           id = '******' then nc need not be an array.
* iyr:   starting year of series
* n:     number of data values in series
*           returned as -999 if end-of-file
* y:     array of index or measurement values
*           when q = 'iw', if maximum value of series is over
*           9.989, entire series is divided by a factor of ten
*           in order to lower the maximum to 9.989 or lower.
* nc:    array of number of series included in index value for each year
*           when q = 'ir', see exception under 'id' above.
*           when q = 'iw', if nc(1) is zero or negative,
*           nc need not be an array)
* ifl:   count of number of series read; returned as -999 if end-of-file
* iu:    unit for reading or writing data
* iuo:   unit for printed output
* q:     character string *2 to indicate action desired:
*
*      indices,      read from file:   q = 'ir'
*      measurements, read from file:   q = 'mr'
*      indices,      write on file:    q = 'iw'
*      measurements, write on file:    q = 'mw'
*
************************************************************
*
      subroutine trrw(id,iyr,n,y,nc,ifl,iu,iuo,q)
      integer iyr,n,ifl,iu,iuo
      character id*(*),ida*8,line*80,q*2
      integer nc(1),iy(10),inc(10),flag
      real y(1)
c      flag=-9999
      if(q .eq. 'ir')then
      k=1
      if(id(1:6) .eq. '******')k=-1
      elseif(q .eq. 'mw' .or. q .eq. 'iw')then
      goto 5
      endif
    2 n=0
      isw=0
      ier=0
      ifl=max0(ifl,0)
      ifl=ifl+1
      goto 1
   91 ier=ier+1
      if(ier .gt. 25)stop'sbr trrw:  over 25 text lines'
c      if(iuo .gt. 0)write(iuo,'(/t10,''title of data: '',a/)')line
c      write(*,'(/''title of data:''/a)')line
    1 read(iu,'(a80)',err=99,end=99)line
      if(n .eq. 0 .and. line(80:80) .eq. '#')then
      backspace iu
      call datar(line,n,y,iyr,iu)
      id=line(1:8)
c      if(iuo .gt. 0)write(iuo,'(/t10,''title of data:  '',a/)')line
      return
      endif
      if(q .eq. 'ir')then
      read(line,'(a,t7,i4,10(i4,i3))',err=91,end=99)
     +ida,iyr1,(iy(i),inc(i),i=1,10)
      ida(7:8)='  '
      elseif(q .eq. 'mr')then
      read(line,'(a,t9,i4,10i6)',err=91,end=99)ida,iyr1,iy
      endif
      isw=isw+1
      l=10
      if(isw .eq. 1)then
      id=ida
      iyr=iyr1
      if(iyr .ge. 0)then
      l=10-mod(iyr,10)
      else
      l=abs(mod(iyr,10))
      if(l .eq. 0)l=10
      endif
      endif
      if(q .eq. 'mr')goto 3
      l=11-l
      do 10 i=l,10
      if(iy(i) .lt. 9990)then
      n=n+1
      y(n)=float(iy(i))*.001
      if(k .ge. 0)nc(n)=inc(i)
      else
      return
      endif
   10 continue
      goto 1
   99 n=-999
      ifl=-999
      return
    3 do 20 i=1,l
      if(iy(i) .eq. -9999)then
        flag=-9999
        do 6 j=1,n
          y(j)=y(j)*.1
    6   continue
        return
      endif
      if(iy(i) .eq. 999)then
        flag=-9999
        return
      endif
      n=n+1
      y(n)=float(iy(i))*.01
   20 continue
      goto 1
    5 if(n .le. 0)return
      ier=0
      j=0
      kyr=iyr
      if(q .eq. 'iw')then
*
* adjust series if max value is over 9.989
      f=1.
      yx=y(1)
      do 60 i=2,n
   60 yx=max(y(i),yx)
   61 if((yx*f) .gt. 9.98949)then
      f=f*.1
      goto 61
      endif
      if(f .ne. 1.)then
      write(*,'(''>> sbr trrw: series '',a8,'' max value'',f12.3,
     +'' multiplied by'',f10.6)')id,yx,f
      endif
      if(iyr .ge. 0)then
      l=mod(iyr,10)+1
      else
      l=11-abs(mod(iyr,10))
      if(l .eq. 11)l=1
      endif
   32 do 30 i=1,10
      iy(i)=9990
   30 inc(i)=0
      do 40 i=l,10
      j=j+1
      if(j .gt. n)then
      j=-999
      goto 41
      endif
      if(nc(1) .gt. 0)then
      inc(i)=nc(j)
      else
      inc(i)=1
      endif
   40 iy(i)=nint(y(j)*1000.*f)
   41 write(iu,'(a6,i4,10(i4,i3))')
     +id(1:6),kyr,(iy(k),inc(k),k=1,10)
      if(j .lt. 0)return
      kyr=kyr-l+11
      l=1
      goto 32
      elseif(q .eq. 'mw')then
      goto 4
      endif
    4 if(iyr .ge. 0)then
      l=10-mod(iyr,10)
      else
      l=abs(mod(iyr,10))
      if(l .eq. 0)l=10
      endif
   42 do 50 i=1,l
      j=j+1
      if(j .gt. n)then
ccc      iy(i)=999
        iy(i)=flag
        l=i
        j=-999
        goto 51
      else
        if(flag.eq.999)then
          iy(i)=nint(y(j)*100.)
          if(iy(i) .eq. 999)iy(i)=1000
        elseif(flag.eq.-9999)then
          iy(i)=nint(y(j)*1000.)
        endif      
      endif
   50 continue
   51 write(iu,'(a8,i4,10i6)')id,kyr,(iy(k),k=1,l)
      if(j .lt. 0)return
      kyr=kyr+l
      l=10
      goto 42
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr datar: read data series preserving significant digits
*
*  programmed may 1986 by
*
*     richard  l  holmes
*     tree-ring laboratory
*     university of arizona
*     tucson, arizona  85721  usa
*
*  modified by r l holmes 19 sep 1988
*
*  call dataw (id,n,x,iyr,iu) write data series
*  call datar (id,n,x,iyr,iu) read data series
*
* parameters:
*
* id:  character string of flexible length to identify series
* n:   number of data values
*        (note: in calling datar, if 'n' is returned as a negative
*         value, an end-of-file has been reached.)
* x:   array of data values
* iyr: year of first data value
* iu:  unit to write or read series
*
      subroutine datar(id,n,x,iyr,iu)
      integer n,iyr,iu
      character line*80,frm*8,s*2,id*(*)  !emt*8,
      real x(1)
      s='rd'
      lz=min(len(id),80)
    2 read(iu,'(a)',err=99,end=99)line
      if(line(1:3) .ne.'n  ')then
      id(1:lz)=line(1:lz)
      endif
      read(line,'(2(2x,i8),t23,a)',err=2,end=99)n,iyr,frm
      if(index(frm,'(') .eq. 0 .or. index(frm,')') .eq. 0)goto 2
      read(iu,frm,err=97,end=98)(x(i),i=1,n)
      return
   97 write(*,'(''>> sbr datar:  error reading data'')')
      return
   98 write(*,'(''>> sbr datar:  unexpected end reading data'')')
      return
   99 n=-999
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr comint: determines optimum common interval
*
*  programmed apr 1986 by
*
*     richard  l  holmes
*     tree-ring laboratory
*     university of arizona
*     tucson, arizona  85721  usa
*
*  modified by r l holmes 22 aug 1988
*  modified by pjk 2/2006
*
      subroutine comint(mxs,iyr,id,n,ns,ou)
      integer mxs,ns,ou
      integer n(mxs),iyr(mxs)
      character*8 id(mxs)
c
c local decs
	integer mxs2,i,l,m,j,k,ifl,icy,lcy,nly
	integer lyr,iz,lly,is,ici,nyc,ny,nci,nsci,nrci
	parameter (mxs2=50000)   !must = mxrec
       integer l20(mxs2,20)
	character sqc(mxs2)*40,sa*40,sb*40
c
*
      if(ns.gt.mxs2)goto 100
      ifl=0
c
      icy= 32000
      lcy=-32000
*
* here is where you would dump your 
* file data into the character array sqc 
      do 10 m=1,ns
          lyr=iyr(m)+n(m)-1
          icy=min(iyr(m),icy)
          lcy=max(lyr,lcy)
          iz=9000-iyr(m)
c          write(sqc(ns),'(5i5,2x,a8,5x)')lyr,iz,n(m),iyr(m),ns,id(m)
          write(sqc(m),'(5i5,2x,a8,5x)')lyr,iz,n(m),iyr(m),ns,id(m)
c          write(*,'(5i5,2x,a8,5x)')lyr,iz,n(m),iyr(m),ns,id(m)
c          pause '5'
10 	continue
c
       do 20 i=1,ns
          sa=sqc(i)
          sb=sqc(i)
          k=i
          do 30 j=i,ns
             if(sqc(j) .gt. sa)then
                sa=sqc(j)
                k=j
             endif
   30     continue
          sqc(k)=sb
20     sqc(i)=sa
*
* make array of starting and ending years of series
       nly=1
c            pause'1-commint'
       read(sqc(1),'(i5,t1,i5,t16,i5)')l20(1,3),l20(1,2),l20(1,1)
       lly=l20(nly,3)
       do 60 i=2,ns
          read(sqc(i),'(i5,t1,i5,t16,i5)',err=1000)l20(mxs,3),l20(i,2),l20(i,1)
          if(lly .ne. l20(mxs,3))then
             nly=nly+1
             l20(nly,3)=l20(mxs,3)
             lly=l20(nly,3)
          endif
60     continue
c	write(*,*)nly
c	pause '2-commint'
       nrci=0
       do 40 l=1,nly
          do 50 i=1,ns
             nc=0
             ny=l20(l,3)-l20(i,1)+1
             do 70 is=1,ns
                 if(l20(is,1) .le. l20(i,1) .and. l20(is,2) .ge. l20(l,3))then
                    nc=nc+1
                 endif
70           continue
             nyc=nc*ny
             if(nyc .gt. nrci)then
                 nrci=nyc
                 ici=l20(i,1)
                 lci=l20(l,3)
                 nci=ny
                 nsci=nc
             endif
50        continue
c	   write(*,*)l,nly,i,ns
c	   pause '3-commint'
40     continue
c	pause '4-commint'
*
* print results on screen
       write(*,'(/3x,42(''*''))')
       write(*,'(2x,'' **  Time span of chronology: '',i5,'' to'',i5)')
     +        icy,lcy
       write(*,'(2x,'' **  Number of series in chronology: '',i6)')ns
       write(*,'(2x,'' **  Optimum common interval: '',i5,'' to'',i5)')
     +        ici,lci
       write(*,'(2x,'' **  Number of years included: '',i12)')nci
       write(*,'(2x,'' **  Number of Series included:     '',i7,/,
     + 2x,'' **  Number of rings:               '',i7)')nsci,nrci
       write(*,'(3x,42(''*''))')
       return
100    write(*,'('' >> file contains more than program limit of'',i4,
     +'' series'')')mxs2-1
       stop' too many series for subroutine Commint'
1000   write(*,*)sqc(i)
	write(*,'(1x,''error reading sqc)'')')
	return
       end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine supsmu (n,x,y,w,iper,span,alpha,smo,sc)
c
c------------------------------------------------------------------
c
c super-smoother.
c
c Friedman J.H. (1984). A variable span smoother. Department of Statistics,
c    Stanford University, Technical Report LCS5.
c
c version 10/10/84.
c
c coded  and copyright (c) 1984 by:
c
c                        Jerome H. Friedman
c                     Department of Statistics
c                               and
c                Stanford Linear Accelerator Center
c                        Stanford University
c
c all rights reserved.
c
c
c input:
c    n : number of observations (x,y - pairs).
c    x(n) : ordered abscissa values.
c    y(n) : corresponding ordinate (response) values.
c    w(n) : weight for each (x,y) observation.
c    iper : periodic variable flag.
c       iper=1 => x is ordered interval variable.
c       iper=2 => x is a periodic variable with values
c                 in the range (0.0,1.0) and peroid 1.0.
c    span : smoother span (fraction of observations in window).
c           span=0.0 => automatic (variable) span selection.
c    alpha : controles high frequency (small span) penality
c            used with automatic span selection (bass tone control).
c            (alpha.le.0.0 or alpha.gt.10.0 => no effect.)
c output:
c   smo(n) : smoothed ordinate (response) values.
c scratch:
c   sc(n,7) : internal working storage.
c
c note:
c    for small samples (n < 40) or if there are substantial serial
c    correlations between obserations close in x - value, then
c    a prespecified fixed span smoother (span > 0) should be
c    used. reasonable span values are 0.2 to 0.4.
c
c------------------------------------------------------------------
c
      real*4 x(n),y(n),w(n),smo(n),sc(n,7)
      common /spans/ spans(3) /consts/ big,sml,eps
      if (x(n).gt.x(1)) go to 30
      sy=0.0
      sw=sy
      do 10 j=1,n
      sy=sy+w(j)*y(j)
      sw=sw+w(j)
 10   continue
      a=0.0
      if (sw.gt.0.0) a=sy/sw
      do 20 j=1,n
      smo(j)=a
 20   continue
      return
 30   i=n/4
      j=3*i
      scale=x(j)-x(i)
 40   if (scale.gt.0.0) go to 50
      if (j.lt.n) j=j+1
      if (i.gt.1) i=i-1
      scale=x(j)-x(i)
      go to 40
 50   vsmlsq=(eps*scale)**2
      jper=iper
      if (iper.eq.2.and.(x(1).lt.0.0.or.x(n).gt.1.0)) jper=1
      if (jper.lt.1.or.jper.gt.2) jper=1
      if (span.le.0.0) go to 60
      call smooth (n,x,y,w,span,jper,vsmlsq,smo,sc)
      return
 60   do 70 i=1,3
      call smooth (n,x,y,w,spans(i),jper,vsmlsq,sc(1,2*i-1),sc(1,7))
      call smooth (n,x,sc(1,7),w,spans(2),-jper,vsmlsq,sc(1,2*i),h)
 70   continue
      do 90 j=1,n
      resmin=big
      do 80 i=1,3
      if (sc(j,2*i).ge.resmin) go to 80
      resmin=sc(j,2*i)
      sc(j,7)=spans(i)
 80   continue
      if (alpha.gt.0.0.and.alpha.le.10.0.and.resmin.lt.sc(j,6).and.resmi
     1n.gt.0.0) sc(j,7)=sc(j,7)+(spans(3)-sc(j,7))*amax1(sml,resmin/sc(j
     2,6))**(10.0-alpha)
 90   continue
      call smooth (n,x,sc(1,7),w,spans(2),-jper,vsmlsq,sc(1,2),h)
      do 110 j=1,n
      if (sc(j,2).le.spans(1)) sc(j,2)=spans(1)
      if (sc(j,2).ge.spans(3)) sc(j,2)=spans(3)
      f=sc(j,2)-spans(2)
      if (f.ge.0.0) go to 100
      f=-f/(spans(2)-spans(1))
      sc(j,4)=(1.0-f)*sc(j,3)+f*sc(j,1)
      go to 110
 100  f=f/(spans(3)-spans(2))
      sc(j,4)=(1.0-f)*sc(j,3)+f*sc(j,5)
 110  continue
      call smooth (n,x,sc(1,4),w,spans(1),-jper,vsmlsq,smo,h)
      return
      end
      subroutine smooth (n,x,y,w,span,iper,vsmlsq,smo,acvr)
      dimension x(n),y(n),w(n),smo(n),acvr(n)
      integer in,out
      double precision wt,fbo,fbw,xm,ym,tmp,var,cvar,a,h,sy
      xm=0.0
      ym=xm
      var=ym
      cvar=var
      fbw=cvar
      jper=iabs(iper)
      ibw=0.5*span*n+0.5
      if (ibw.lt.2) ibw=2
      it=2*ibw+1
      do 20 i=1,it
      j=i
      if (jper.eq.2) j=i-ibw-1
      xti=x(j)
      if (j.ge.1) go to 10
      j=n+j
      xti=x(j)-1.0
 10   wt=w(j)
      fbo=fbw
      fbw=fbw+wt
      if (fbw.gt.0.0) xm=(fbo*xm+wt*xti)/fbw
      if (fbw.gt.0.0) ym=(fbo*ym+wt*y(j))/fbw
      tmp=0.0
      if (fbo.gt.0.0) tmp=fbw*wt*(xti-xm)/fbo
      var=var+tmp*(xti-xm)
      cvar=cvar+tmp*(y(j)-ym)
 20   continue
      do 80 j=1,n
      out=j-ibw-1
      in=j+ibw
      if ((jper.ne.2).and.(out.lt.1.or.in.gt.n)) go to 60
      if (out.ge.1) go to 30
      out=n+out
      xto=x(out)-1.0
      xti=x(in)
      go to 50
 30   if (in.le.n) go to 40
      in=in-n
      xti=x(in)+1.0
      xto=x(out)
      go to 50
 40   xto=x(out)
      xti=x(in)
 50   wt=w(out)
      fbo=fbw
      fbw=fbw-wt
      tmp=0.0
      if (fbw.gt.0.0) tmp=fbo*wt*(xto-xm)/fbw
      var=var-tmp*(xto-xm)
      cvar=cvar-tmp*(y(out)-ym)
      if (fbw.gt.0.0) xm=(fbo*xm-wt*xto)/fbw
      if (fbw.gt.0.0) ym=(fbo*ym-wt*y(out))/fbw
      wt=w(in)
      fbo=fbw
      fbw=fbw+wt
      if (fbw.gt.0.0) xm=(fbo*xm+wt*xti)/fbw
      if (fbw.gt.0.0) ym=(fbo*ym+wt*y(in))/fbw
      tmp=0.0
      if (fbo.gt.0.0) tmp=fbw*wt*(xti-xm)/fbo
      var=var+tmp*(xti-xm)
      cvar=cvar+tmp*(y(in)-ym)
 60   a=0.0
      if (var.gt.vsmlsq) a=cvar/var
      smo(j)=a*(x(j)-xm)+ym
      if (iper.le.0) go to 80
      h=0.0
      if (fbw.gt.0.0) h=1.0/fbw
      if (var.gt.vsmlsq) h=h+(x(j)-xm)**2/var
      acvr(j)=0.0
      a=1.0-w(j)*h
      if (a.le.0.0) go to 70
      acvr(j)=abs(y(j)-smo(j))/a
      go to 80
 70   if (j.le.1) go to 80
      acvr(j)=acvr(j-1)
 80   continue
      j=1
 90   j0=j
      sy=smo(j)*w(j)
      fbw=w(j)
      if (j.ge.n) go to 110
 100  if (x(j+1).gt.x(j)) go to 110
      j=j+1
      sy=sy+w(j)*smo(j)
      fbw=fbw+w(j)
      if (j.lt.n) go to 100
 110  if (j.le.j0) go to 130
      a=0.0
      if (fbw.gt.0.0) a=sy/fbw
      do 120 i=j0,j
      smo(i)=a
 120  continue
 130  j=j+1
      if (j.le.n) go to 90
      return
      end
      block data
      common /spans/ spans(3) /consts/ big,sml,eps
c
c---------------------------------------------------------------
c
c this sets the compile time (default) values for various
c internal parameters :
c
c spans : span values for the three running linear smoothers.
c spans(1) : tweeter span.
c spans(2) : midrange span.
c spans(3) : woofer span.
c (these span values should be changed only with care.)
c big : a large representable floating point number.
c sml : a small number. should be set so that (sml)**(10.0) does
c       not cause floating point underflow.
c eps : used to numerically stabilize slope calculations for
c       running linear fits.
c
c these parameter values can be changed by declaring the
c relevant labeled common in the main program and resetting
c them with executable statements.
c
c-----------------------------------------------------------------
c
      data spans,big,sml,eps /0.05,0.2,0.5,1.0e20,1.0e-7,1.0e-3/
      end
*********************************************************************************
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c function to find the maximum value
	function maximum(iarray, nsize)
	dimension iarray(nsize)
	maximum=iarray(1)
	do 10 j=2,nsize
	   if (maximum.lt.iarray(j)) maximum=iarray(j)
10     continue
	return
	end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr sort
*
*      heapsort algorithm from numerical recipes
*
*      n   =  number of observations in ra to be sorted
*      ra  =  input series to be sorted. ra is returned
*             with its observations ranked in ascending
*             numerical order.
*
      subroutine sort(n,ra)
      real*4 ra(n)
      l=n/2+1
      ir=n
   10 continue
        if(l.gt.1)then
          l=l-1
          rra=ra(l)
        else
          rra=ra(ir)
          ra(ir)=ra(1)
          ir=ir-1
          if(ir.eq.1)then
            ra(1)=rra
            return
          endif
        endif
        i=l
        j=l+l
   20   if(j.le.ir)then
          if(j.lt.ir)then
            if(ra(j).lt.ra(j+1))j=j+1
          endif
          if(rra.lt.ra(j))then
            ra(i)=ra(j)
            i=j
            j=j+j
          else
            j=ir+1
          endif
        go to 20
        endif
        ra(i)=rra
      go to 10
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

