	program SignalFree
c
	integer iu,ou,mxd,mxr
	parameter(mxd=10000,mxr=10000,iu=69,ou=68)
	integer p,j,k,ifl,ls,nrec
	integer mniyr,mxlyr,mnage,mxage,ici,lci,str(mxr),stp(mxr)
	integer age(mxr),iyr(mxr),lyr(mxr),numsft(mxr)
	real*4   rrw(mxr,mxd),sft(mxr,mxd)
	real*4   y1(mxd),y2(mxd),nc(mxd),nrrw(mxd),sftindx(mxr,mxd)
	real*4   crn1(mxd),crn2(mxd),sum(mxd),sumsft,sumrrw
	real*4   ones(mxd)        !function the returns the sum of array elements
	character*1 z
	character*8 id(mxr)
	character*40 infl,outfl
c
	write(*,'(//1x,''Enter input file name : '',$)')
	read(*,'(a)')infl
	do p=1,40-1
	  read(infl(p:p),'(a)')z
	  if(z.eq.' ')goto 5
	end do
5	outfl=infl(1:p-1)//'_out'
	write(*,'(/1x,''input file  = '',a)')infl
	write(*,'(1x,''output file = '',a)')outfl
	write(*,'(//1x,''Compute common period'')')
	write(*,'(1x,'' [f]ancypants'')')
	write(*,'(1x,'' [s]imple     : '',$)')
	read(*,'(a)')z
c
	mniyr=mxd
	mxlyr=-1.0*mniyr
	ls=100
c open and read
	open(unit=iu,file=infl,status='old',readonly)
	p=1
10	call trrw(mxd,id(p),iyr(p),age(p),y1,nc,ifl,iu,iuo,'MR') 
c	  write(*,'(1x,''ifl= '',i5)')ifl
	  if(ifl.gt.0)then 
	    nrec=ifl
           lyr(p)=iyr(p)+age(p)-1
           k=1
	    do j=iyr(p),lyr(p)
	      rrw(p,j)=y1(k)
	      k=k+1
	    end do	    
	    if(iyr(p).lt.mniyr)mniyr=iyr(p)
	    if(lyr(p).gt.mxlyr)mxlyr=lyr(p)
	    if(age(p).lt.mnage)mnage=age(p)
	    if(age(p).gt.mxage)mxage=age(p)
	    p=p+1
	    goto 10
	  else
	    write(*,'(/1x,''INPUT file read without errors'')')
	  endif
c
	write(*,'(/1x,'' Number of series read = '',i5)')nrec
	write(*,'(1x,''series   id        iyr   lyr   age   '')')
	do p=1,nrec
	   write(*,'(1x,i5,3x,a8,3i6)')p,id(p),iyr(p),lyr(p),age(p)
 	end do
 	close(iu)
c find common period  this may not be the best interval, maybe a better is one where n.ge.2
	if(z.eq.'F'.or.z.eq.'f')then
	  call comint(mxr,iyr,id,age,nrec,ici,lci)
	else
	  do p=1,nrec
	    do k=iyr(p),lyr(p)
	      nc(k)=nc(k)+1.0
	    end do
	  end do
	  do k=mniyr,mxlyr
	    if(nc(k).ge.2.0)then
	      ici=k
	      goto 11
	    end if
	  end do
11	  do k=mxlyr,mniyr,-1
	    if(nc(k).ge.2.0)then
	      lci=k
	      goto 12
	    endif
	  end do
12	  write(*,'(1x,2i6)')ici,lci 
	endif
c determine where start and stop indices for each series arrays are within common period
	do p=1,nrec
	   do k=ici,lci
	     if(iyr(p).le.k)then
	       str(p)=k
	       goto 20
	     end if
	   end do
20	   do k=lci,ici,-1
	     if(lyr(p).ge.k)then
	       stp(p)=k
	       goto 30
	     end if
	   end do
30	   numsft(p)=stp(p)-str(p)+1
	end do
c
	do p=1,mxd
	  sum(p)=0.00
	  nrrw(p)=0.00
	  ones(p)=1.0
	end do
c
c spline all rrw series initial chronolonogy
	do p=1,nrec
	   k=1
	   do j=iyr(p),lyr(p)
	      y1(k)=rrw(p,j)
	      k=k+1
	   end do
	   call spline(mxd,age(p),y1,y2,ls)
	   k=1
	   do j=iyr(p),lyr(p)
	      sftindx(p,j)=y1(k)/y2(k)
	      sum(j)=sum(j)+sftindx(p,j)  !sum indices
	      nrrw(j)=nrrw(j)+1
	      k=k+1
	   end do
	end do
	do p=mniyr,mxlyr
	   crn1(p)=sum(p)/nrrw(p)  !original chronology total length of longest series
	   crn2(p)=crn1(p)         ! a copy
	end do
c
40	continue   !begin loop until convergence of some kind or another
c
	  do k=mniyr,mxlyr
	    sum(k)=0.00
	  end do
	  do p=1,nrec  !for all trees 
	    write(*,*)str(p),stp(p),ici,lci
	    do k=str(p),stp(p)         
	        sft(p,k)=rrw(p,k)/crn2(k)   ! 2)remove currently estimated common signal where n.ge.2
	        y1(k)=rrw(p,k)
	        y2(k)=sft(p,k)
	    end do
	    sumrrw=asum(mx,y1,str(p),stp(p)) !sum all rrw for tree p
	    sumsft=asum(mx,y2,str(p),stp(p)) !sum all sftree's for tree p
	    write(*,'(1x,2f10.4)')sumrrw,sumsft
	    do k=str(p),stp(p)
	       sft(p,k)=sft(p,k)*(sumrrw/sumsft)  ! 3) reset means
	    end do
	    sumrrw=0.00
	    sumsft=0.00
	    j=1
	    do k=str(p),stp(p)
	       y1(j)=sft(p,k)
	       j=j+1
	    end do
	    call spline(mxd,numsft(p),y1,y2,ls)
	    j=1
	    do k=str(p),stp(p)
	         sftindx(p,k)=y1(j)/y2(j)  !  4) standardize all sft's
	         sum(k)=sum(k)+sftindx(p,k)
	         j=j+1
	    end do	    
	  end do
c make new corrected signal free crn
	  do k=ici,lci
	    crn2(k)=crn2(k)*(sum(k)/nrrw(k))
	  end do
c
c plot something
	  do k=ici,lci
	    write(*,'(1x,i4,3f8.3)')k,crn1(k),crn2(k)
	  end do
	  call display(mxd,mxr,nrec,rrw,crn1,crn2,ones,ici,lci)

c go back to top and continue removing  signal
	  write(*,'(/1x,''Do you want to continue? '',$)')
	  write(*,'(1x,''[Y]es or [N]o '',$)')
	  read(*,'(a)')z
	  if(z.eq.'Y'.or.z.eq.'y'.or.z.eq.' ')then
	    goto 40
	  else
	    goto 100
	  endif
	go to 40
c
100	write(*,'(/1x,''Normal termination of program'')')
	stop
	end
c**********************************************************
* sum input array y from str to stp: 
*
       function asum(mx,y,istr,istp)
       integer mx,istr,istp,p
       real y(mx)
       asum=0.0
	do p=istr,istp
	  asum=asum+y(p)
	end do
	return
	end
c****
c**************************************************************
c Subroutine display, plots bai computed  estimates 
c pjk-trl, ldeo 3/2007
c
      subroutine display(mxd,mxr,nrec,rrw,crn1,crn2,crn3,cyr1,cyr2)
c
      integer mxd,mxr,cyr1,cyr2,nrec
	real*4   rrw(mxr,mxd)
	real*4   crn1(mxd),crn3(mxd),crn2(mxd)
c
c local declarations
      integer   mx
      parameter (mx=10000)    !local parameters must equal mxd and mxr in main program
      integer   p,j,k,n1
      integer*4 ifault,mnbtic
      real*4    zhigh,zlow,y3high,y3low,mxbtic
      real*4    y1(mx),y2(mx),y3(mx),x1(mx)
      character*2  ans
c
c *** plotstyle variables
	integer*4 gopt,opt,mxvar
	parameter (mxvar=20)
	integer*4 red,green,blue,pen(mxvar),colline(mxvar),coltext,colbox
	integer*4 mnxtic,mnytic,mnztic,symbol(mxvar),colsym(mxvar)
	real*4    mxxtic,mxytic,mxztic,xhigh,xlow,y1high,y1low,y2high,y2low
       character*8 dev
c      
c begin routine: set device: dev='Mac1'(terminal), dev='psc'(color postscript), dev='ps'(monochrome ps)
c                dev='win3' for pc displays
c
10     gopt=99
       call pltstyle(dev,gopt,opt,mxvar,red,green,blue,pen,colline,
     *       coltext,colbox,mnxtic,mnytic,mnztic,symbol,colsym,
     *       mxxtic,mxytic,mxztic,xhigh,xlow,y1high,y1low,
     *       y2high,y2low) 
c
       dev='Mac1'
c
	n1=cyr2-cyr1+1
	k=1
       do j=cyr1,cyr2
             x1(k)=real(j)
             y1(k)=crn1(j)
             y2(k)=crn2(j)
             y3(k)=crn3(j)
             if(x1(k).gt.xhigh)xhigh=x1(k)
             if(x1(k).lt.xlow)xlow=x1(k)
             k=k+1
       end do
c find x and y scale axis parameters
	call SCALE(xlow, xhigh, 10, zlow, mxxtic, zhigh, IFAULT)
	xlow=zlow;xhigh=zhigh
	mnxtic=10
	y1low=-1.0
	y1high=3.0
	mnytic=2
	mxytic=1.0
c	    
c set device: dev='Mac1'(terminal), dev='psc'(color postscript), dev='ps'(monochrome ps)
35	call plsdev(dev)
c set background color
	    call plscol0(0,255,255,255)              !change default color map black>white 
	    call plscol0(15,0,0,0)                   !change default color map white>black
	    call plscolbg(red,green,blue)	      !change background color if other than white is desired
c    initialize plplot
	    call plinit()
	    call plfontld(1)
	    call plfont(1)
c    define the viewport one column by one row
	    call plssub(1,1)
	    call pladv(0)
	    call plvpor(0.12,0.90,0.2,0.85)   !plot size
	    call plwind(xlow,xhigh,y1low,y1high)
c set box color
	    call plcol(colbox)   !color box and lables
	    call plaxes (xlow,y1low,'abinst',mxxtic,mnxtic,'abinstv',mxytic,mnytic)
c set lable color
	    call plcol(coltext)   !color box and lables
	    call pllab('YEAR','Index','SignalFree')  !plot raw ring width
c plot lines
c set line style
	    call pllsty(pen(1))
c set line width
           call plwid(2)
c set line color
	    call plcol(colline(1))  !color lines (black)
c plot series
           call plline(n1,x1,y1)
	    call plmtex('b',-3.2,0.14,0.0,' Original')
c set line width
           call plwid(1)
c set line color
	    call plcol(colline(2))  !color lines (red)
c set line style
	    call pllsty(pen(2))
c plot series
           call plline(n1,x1,y2)
	    call plmtex('b',-3.2,0.54,0.0,'Corrected Original')
c set line color
	    call plcol(colline(3))  !color lines (red)
c set line style
	    call pllsty(pen(1))
c plot series
           call plline(n1,x1,y3)
c	    call plmtex('b',-3.2,0.74,0.0,'blank')
c
c end plotting 
	call plend()
      opt=1
      ans='c '
      dev='Mac1'
51    write(*,'(//15x,''|************** plot options menu **************|'')')
      write(*,'(15x,''| replot             [0] max. x            [14] |'')')
      write(*,'(15x,''| background color   [1] min. x            [15] |'')')
      write(*,'(15x,''| frame color        [2] max. y1           [16] |'')')
      write(*,'(15x,''| text color         [3] min. y1           [17] |'')')
      write(*,'(15x,''| major ticks x      [8]                        |'')')
      write(*,'(15x,''| minor ticks x      [9]                        |'')')
      write(*,'(15x,''| major ticks y1    [10]                        |'')')
      write(*,'(15x,''| minor ticks y1    [11]                        |'')')
      write(*,'(15x,''|                                               |'')')                           
      write(*,'(15x,''|                                               |'')')                         
      write(*,'(15x,''|                                               |'')')
      write(*,'(15x,''|                reset options  [99]            |'')')
      write(*,'(15x,''|-----------------------------------------------|'')')
      write(*,'(15x,''|                [c]ontinue  [s]top             |'')') 
      write(*,'(15x,''|_______________________________________________|'')')
      write(*,'(/15x,'' enter option number > '',$)')
      read(*,'(a2)')ans
      if((ans.eq.'C ').or.(ans.eq.'c ').or.(ans.eq.'  ').or.
     *          (ans.eq.'Y ').or.(ans.eq.'y '))then
         return
      elseif((ans.eq.'S ').or.(ans.eq.'s '))then
         stop
      else
         read(ans(1:2),'(i2)')gopt     !convert ans to gopt then continue
      endif
c
      if(gopt.eq.99)goto 10
      if(gopt.eq.0)then
         goto 35
      elseif((gopt.gt.0).and.(gopt.le.19))then
         call pltstyle(dev,gopt,opt,mxvar,red,green,blue,pen,colline,
     *       coltext,colbox,mnxtic,mnytic,mnztic,symbol,colsym,
     *       mxxtic,mxytic,mxztic,xhigh,xlow,y1high,y1low,
     *       y2high,y2low)                          
         goto 51
      else
         write(*,'(/15x,''Unknown command - try again'')')
         pause'               Hit return to continue'
         goto 51
      endif
      return
cc
      end
c***********************************************************************
      SUBROUTINE SCALE(FMN, FMX, N, VALMIN, STEP, VALMAX, IFAULT)
C
C     ALGORITHM AS 96  APPL. STATIST. (1976) VOL.25, NO.1
C
C     Given extreme values FMN, FMX, and the need for a scale with N
C     marks, calculates value for the lowest scale mark (VALMIN) and
C     step length (STEP) and highest scale mark (VALMAX).
C
      REAL FMN, FMX, VALMIN, STEP, VALMAX
      INTEGER N, IFAULT
C
C     Units for step lengths
C
      REAL UNIT(11)
C
C     Local variables
C
      INTEGER NUNIT, I, J
      REAL TOL, ZERO, HALF, ONE, TEN, BIAS, FMAX, FMIN, RN, X, S, RANGE
C
C     Array length unit()
C
      DATA NUNIT /11/
C
C     Local constant, defining effective equality of values.
C
      DATA TOL /5.0E-6/
      DATA ZERO /0.0/, HALF /0.5/, ONE /1.0/, TEN /10.0/
      DATA BIAS /1.0E-4/
      DATA UNIT /1.0, 1.2, 1.6, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0,
     *        10.0/
C
      FMAX = FMX
      FMIN = FMN
      IFAULT = 1
C
C     Test for valid parameter values
C
      IF (FMAX .LT. FMIN .OR. N .LE. 1) RETURN
      IFAULT = 0
      RN = N - 1
      X = ABS(FMAX)
      IF (X .EQ. ZERO) X = ONE
      IF ((FMAX - FMIN) / X .GT. TOL) GO TO 20
C
C     All values effectively equal
C
      IF (FMAX .LT. ZERO) THEN
        FMAX = ZERO
      ELSE IF (FMAX .EQ. ZERO) THEN
        FMAX = ONE
      ELSE
        FMIN = ZERO
      END IF
C
   20 STEP = (FMAX - FMIN) / RN
      S = STEP
C
C     Find power of 10
C
   25 IF (S .GE. ONE) GO TO 30
      S = S * TEN
      GO TO 25
   30 IF (S .LT. TEN) GO TO 35
      S = S / TEN
      GO TO 30
C
C     Calculate STEP
C
   35 X = S - BIAS
      DO 40 I = 1, NUNIT
        IF (X .LE. UNIT(I)) GO TO 45
   40 CONTINUE
   45 STEP = STEP * UNIT(I) / S
      RANGE = STEP * RN
C
C     Make first estimate of VALMIN
C
      X = HALF * (ONE + (FMIN + FMAX - RANGE) / STEP)
      J = X - BIAS
      IF (X .LT. ZERO) J = J - 1
      VALMIN = STEP * FLOAT(J)
C
C     Test if VALMIN could be zero
C
      IF (FMIN .GE. ZERO .AND. RANGE .GE. FMAX) VALMIN = ZERO
      VALMAX = VALMIN + RANGE
C
C     Test if VALMAX could be zero
C
      IF (FMAX .GT. ZERO .OR. RANGE .LT. -FMIN) RETURN
      VALMAX = ZERO
      VALMIN = -RANGE
      RETURN
      END
cSUBROUTINE PLOTSTYLE *******************************************************************
	subroutine pltstyle(dev,gopt,opt,mxvar,red,green,blue,pen,colline,
     *       coltext,colbox,mnxtic,mnytic,mnztic,symbol,colsym,
     *       mxxtic,mxytic,mxztic,xhigh,xlow,yhigh,ylow,
     *       zhigh,zlow)             
c
c	sets up variables for passing to plplot routines that govern color and 
c	line style changes to plots and all sorts of bells and whistles -pjk
c
	integer*4 gopt,opt,mxvar
	integer*4 red,green,blue,pen(mxvar),colline(mxvar),coltext,colbox
	integer*4 mnxtic,mnytic,mnztic,symbol(mxvar),colsym(mxvar)
	real*4 mxxtic,mxytic,mxztic,xhigh,xlow,yhigh,ylow,zhigh,zlow
      character*8 dev
c local
c	integer*4 winu2
c	parameter(winu2=92)
	integer*4 tmp
c
c Begin subroutine
	if((gopt.gt.1.and.gopt.le.5).or.(gopt.eq.21))then
c	  open(unit=*,file='PlotStyleWindow',access='window,290,522,26,388')
	  write(*,'(/13x,''color table ______________________________________________'')')
	  write(*,'(13x,''white(0) red(1)  yellow(2)  green(3)  aquamarine(4) pink(5) '')')
	  write(*,'(13x,''wheat(6)  grey(7)  browm(8)  blue(9)  blueviolet(10)'')')
	  write(*,'(13x,''cyan(11)  turquoise(12)  magenta(13)  salmon(14)  black(15)'')')
	  write(*,'(13x,''___________________________________________________________'',/)')
	else if(gopt.eq.1)then
c	  open(unit=*,file='plotstylewindow',access='window,290,522,26,388')
	  write(*,'(/13x,''color table (background only) ____________________________'')')
	  write(*,'(13x,''black(0) white(1)  user defined (2)'')')
	  write(*,'(13x,''___________________________________________________________'',/)')
	else if(gopt.ge.6.and.gopt.le.21)then
c	  open(unit=*,file='PlotStyleWindow',access='window,290,522,26,388')
	endif  
	select case (gopt)
	   case(1)
	      write(*,'(/13x,''select new background color [0/1/2] : '',$)')
	      read(*,*)tmp
	      if(tmp.eq.0)then
	        red=0
	        blue=0
	        green=0
	      elseif(tmp.eq.1)then
	        red=254
	        green=254
	        blue=254
	      elseif(tmp.eq.2)then
	        write(*,'(13x,''enter red value [0-255] > '',$)')
	        read(*,*)red
	        write(*,'(13x,''enter green value [0-255] > '',$)')
	        read(*,*)green
	        write(*,'(13x,''enter blue value [0-255] > '',$)')
	        read(*,*)blue
	      endif
c	      close(winu2)
	      return
	   case(2)
	      write(*,'(13x,''select new frame color [0-15] : '',$)')
	      read(*,*)colbox
c	      close(winu2)
	      return
	   case(3)
	      write(*,'(/13x,''select new text color: [0-15] : '',$)')
	      read(*,*)coltext
c	      close(winu2)
	      return
	   case(4)
	      write(*,'(13x,''select new line color: [0-15] : '',$)')
	      read(*,'(i4)')colline(1)
c	      close(winu2)
	      return
	   case(5)
	      write(*,'(13x,''select new curve color: [0-15] : '',$)')
	      read(*,'(i4)')colline(2)
c	      close(winu2)
	      return
	   case(6)
	      write(*,'(/13x,''select new line style: [0-8] : '',$)')
	      read(*,'(i4)')pen(1)
c	      close(winu2)
	      return
	   case(7)
	      write(*,'(/13x,''select new curve width: [0-8]  : '',$)')
	      read(*,'(i4)')pen(2)
c	      close(winu2)
	      return
	   case(8)
	      write(*,*)
	      write(*,*)'             change x axis major tics from: ',mxxtic
	      write(*,'(13x,''to > '',$)')
	      read(*,*)mxxtic
c	      close(winu2)
	      return
	   case(9)
	      write(*,'(/13x,''change x axis minor tics from: '',i8)')mnxtic
	      write(*,'(13x,''to > '',$)')
	      read(*,*)mnxtic
c	      close(winu2)
	      return
	   case(10)
	      write(*,*)'             change y axis major tics from: ',mxytic
	      write(*,'(13x,''to > '',$)')
	      read(*,*)mxytic
c	      close(winu2)
	      return
	   case(11)
	      write(*,'(/13x,''change y axis minor tics from: '',i8)')mnytic
	      write(*,'(13x,''to > '',$)')
	      read(*,'(i4)')mnytic
c	      close(winu2)
	      return
	   case(12)
	      write(*,*)'             change y2 axis major tics from: ',mxztic
	      write(*,'(13x,''to > '',$)')
	      read(*,*)mxztic
c	      close(winu2)
	      return
	   case(13)
	      write(*,'(/13x,''change y2 axis minor tics from: '',i8)')mnztic
	      write(*,'(13x,''to > '',$)')
	      read(*,'(i4)')mnztic
c	      close(winu2)
	      return
	   case(14)
	      write(*,*)'             change x axis max value from: ',xhigh
	      write(*,'(13x,''to > '',$)')
	      read(*,*)xhigh
c	      close(winu2)
	      return
	   case(15)
	      write(*,*)'             change x axis min value from: ',xlow
	      write(*,'(13x,''to > '',$)')
	      read(*,*)xlow
c	      close(winu2)
	      return
	   case(16)
	      write(*,*)'             change y axis max value from: ',yhigh
	      write(*,'(13x,''to > '',$)')
	      read(*,*)yhigh
c	      close(winu2)
	      return
	   case(17)
	      write(*,*)'             change y axis min value from: ',ylow
	      write(*,'(13x,''to > '',$)')
	      read(*,*)ylow
c	      close(winu2)
	      return
	   case(18)
	      write(*,*)'             change y2 axis max value from: ',zhigh
	      write(*,'(13x,''to > '',$)')
	      read(*,*)zhigh
c	      close(winu2)
	      return
	   case(19)
	      write(*,*)'             change y2 axis min value from: ',zlow
	      write(*,'(13x,''to > '',$)')
	      read(*,*)zlow
c	      close(winu2)
	      return
	   case(20)
	      write(*,'(/13x,''change symbol [0-126, 32 = no symbol]: '',i4)')symbol(1)
	      write(*,'(13x,''to > '',$)')
	      read(*,*)symbol(1)
c	      close(winu2)
	      return
	   case(21)
	      write(*,'(/13x,''change symbol color from: '',i4)')colsym(1)
	      write(*,'(13x,''to > '',$)')
	      read(*,*)colsym(1)
c	      close(winu2)
	      return
	   case(99)     !reset all options
	      red=254
	      blue=254
	      green=254
	      coltext=15
	      colbox=15
	      pen(1)=1
	      pen(2)=2
	      colsym(1)=9
	      symbol(1)=20
	      colline(1)=15
	      colline(2)=1
	      colline(3)=9
	      dev='Mac1'
 	      xlow=1000000.0;xhigh=-1000000.00
 	      ylow=xlow;yhigh=xhigh
 	      zlow=xlow;zhigh=xhigh
	end select   
100	return
	end
c****************************************************************
*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,ici,lci)
      integer mxs,ns,ici,lci
      integer n(mxs),iyr(mxs)
      character*8 id(mxs)
c
c local decs
	integer mxs2
	parameter (mxs2=10000)   !must = mxs
       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'
       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
       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
40     continue
*
* print results on screen
       write(*,'(/2x,42(''*''))')
       write(*,'(1x,'' **  Time span of chronology: '',i5,'' to'',i5)')
     +        icy,lcy
       write(*,'(1x,'' **  Number of series in chronology: '',i6)')ns
       write(*,'(1x,'' **  Optimum common interval: '',i5,'' to'',i5)')
     +        ici,lci
       write(*,'(1x,'' **  Number of years included: '',i12)')nci
       write(*,'(1x,'' **  Number of Series included:       '',i5,/,
     + 1x,'' **  Number of rings:               '',i7)')       nsci,nrci
       write(*,'(2x,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)'')')
	stop
       end
c
c**************************************************************************
c_________________________________________________________________________
c**********************************************************************
*SBR SPLINE: FITS CUBIC SPLINE CURVE TO TIME SERIES *******************
c**********************************************************************
*
*  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(maxdat,N,Y,F,LS)
      integer maxdat,ls,N
c      IMPLICIT REAL*8 (A-H,O-Z)
      integer mxy
      PARAMETER (MXY=10000)
      real*8 A(MXY,4),C1(4),C2(3)
      real*8 pct,pi,v,arg,rn,d1,d2,sum
      REAL*4 Y(maxdat),F(maxdat)
      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=DBLE(LS)
      ARG=(2.D0*PI)/V
      P=(6.D0*(DCOS(ARG)-1.D0)**2.D0)/(DCOS(ARG)+2.D0)
      P1=((1.D0/(1.D0-PCT)-1.D0)*6.D0*(DCOS(PI*2.D0/V)-1.D0)**2)/
     +     (DCOS(PI*2.D0/V)+2.D0)
      P2=(6.D0*(DCOS(ARG)-1.D0)**2.D0)/((2.D0**0.5D0-1.D0)*
     *		(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
c****************************************************
C___________________________________________________________
	
*  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(maxdat,ID,IYR,N,Y,NC,IFL,IU,IUO,Q)
      CHARACTER ID*(*),IDA*8,LINE*80,Q*2
      INTEGER NC(maxdat),IY(10),INC(10),FLAG
      REAL Y(maxdat)
      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'
      IF(IUO .GT. 0)WRITE(IUO,'(/T10,''Title of data: ''
     +,A/)')LINE
      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)
      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
*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)
      CHARACTER LINE*80,FRM*8,S*2,ID*(*)
      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
