*Program ArstanSigFreeTom
*
      Program ArstanSigFreeTom
c
c *** does the minimum and implements Tom Melvin's signal free method
c
      parameter(mxy=5000,mxs=1000)
      parameter(iu10=10,iu11=11,iu12=12)
      real*4 x(mxy,mxs),y(mxy),f(mxy),sfind(mxy,mxs)
      real*4 crn1(mxy,2),crn2(mxy,2)
      real*8 tol
ccc      real*4 ybar(mxy,3),rnss(mxy),runr(mxy,5),crns(mxy,5)
      integer iopt,irob,iproc
      integer idate(mxs,3),idate1(mxs,3),nss(mxy),nss1(mxy),ls(mxy)
      character filnam*40,fil1*50,itrdb(3)*85,ident(mxs)*8,id*8
ccc      character tit*85,xlab*20,ylab*20,annot1*35,annot2*35
      character irob1,itrun,melvin,reset,switch,itol,cproc,tag*3
c
c *** open data file to process
c
      call blank(filnam)
      write(*,'(/1x,''Tree-ring data file name (/ to exit)'',t50,''==> '',$)')
      read(*,'(a)')filnam
      if(filnam(1:1).eq.'/')stop
      open(iu10,file=filnam,status='old')
c
c *** first read in the data
c
      ncrns=0
      do i=1,mxy
        do j=1,mxs
          x(i,j)=-9.99
          sfind(i,j)=-9.99
        enddo
      enddo
      call trread(x,y,nc,idate,ident,itrdb,iu10,mxy,mxs)
      close(iu10)
      ifmny=idate(nc+1,1)
      ilmxy=idate(nc+1,2)
      maxn=ilmxy-ifmny+1
      write(*,'(/t30,''    NS   IFY   ILY   NYR'')')
      write(*,'(1x,''Total series read in:'',t30,4i6,''  Years total time span'')')
     *  nc,ifmny,ilmxy,maxn
c
c *** get unique character string for output files
c
      ll=len(filnam)
      jj=0
      do j=1,ll
        jj=ll-j+1
        if(filnam(jj:jj).ne.' ')goto 1
      enddo
    1 ll=jj
      do j=1,ll
        jj=j
        if(filnam(j:j).eq.'.')then
          jj=j-1
          goto 2
        endif
      enddo
    2 ll=jj
c
c *** now find the minimum, lower hinge, median, upper hinge, and maximum segment lengths
c
    3 do j=1,nc
        f(j)=idate(j,3)
      enddo
      call mfs(f,nc,xmedseg,hhseg,hlseg,fsdseg,rmnseg,rmxseg)
      mnseg=int(rmnseg)
      mlseg=int(hlseg)
      medseg=int(xmedseg)
      mhseg=int(hhseg)
      mxseg=int(rmxseg)
      nsegmin=0
      nseghl=0
      nsegmed=0
      nseghh=0
      nsegmax=0
      do j=1,nc
        nyrs=idate(j,3)
        if(nyrs.eq.mnseg)nsegmin=nsegmin+1
        if(nyrs.ge.mlseg)nseghl=nseghl+1
        if(nyrs.ge.medseg)nsegmed=nsegmed+1
        if(nyrs.ge.mhseg)nseghh=nseghh+1
        if(nyrs.eq.mxseg)nsegmax=nsegmax+1
      enddo
      write(*,'(/1x,''ORIGINAL DATA SET'',t30,''   MIN   25%   MED   75%   MAX'')')
      write(*,'(1x,''Segment length statistics:'',t30,5i6)')mnseg,mlseg,medseg,mhseg,mxseg
      write(*,'(1x,''Series .ge. segment length:'',t30,5i6)')nsegmin,nseghl,nsegmed,nseghh,nsegmax
c
c *** process all the data or censor by segment length or starting year
c
    4 cproc=' '
      write(*,'(/1x,''Series processing options (/ to exit):''/
     *  1x,''All series processed as read in'',t39,''<1>''/
     *  1x,''Minimum segment length series'',t40,''2''/
     *  1x,''Series on or before starting year'',t40,''3''/
     *  1x,''Series over common period only'',t40,''4'',t50,''==> '',$)')
      read(*,'(a)')cproc
      if(cproc.eq.'/')stop
      if(cproc.eq.' ')cproc='1'
c
c *** censor the data by chosen method and save it in Tucson raw ring width format
c
      if(cproc.ne.'1')then
        if(cproc.eq.'2')then
    5     mnsl=0
          write(*,'(/1x,''Enter the minimum segment length'',t50,''==> '',$)')
          read(*,*,err=5)mnsl
          call cenbyseg(mnsl,x,nc,idate,ident,mxy,mxs)
          icen1=mnsl
          icen2=0
        elseif(cproc.eq.'3')then
    6     mxsy=0
          write(*,'(/1x,''Enter the maximum starting year'',t50,''==> '',$)')
          read(*,*,err=6)mxsy
          call cenbyify(mxsy,x,nc,idate,ident,mxy,mxs)
          icen1=mxsy
          icen2=0
        elseif(cproc.eq.'4')then
    7     ifc=0
          ilc=0
          write(*,'(/1x,''Enter the common start year'',t50,''==> '',$)')
          read(*,*,err=7)ifc
          write(*,'(1x,''Enter the common end year'',t50,''==> '',$)')
          read(*,*,err=7)ilc
          call cenbycom(ifc,ilc,x,nc,idate,ident,mxy,mxs)
          icen1=ifc
          icen2=ilc
        else
          goto 4
        endif
        call savcen(cproc,icen1,icen2,x,y,nc,idate,ident,itrdb,filnam,iu11,mxy,mxs)
      endif
c
c *** now find the minimum, lower hinge, median, upper hinge, and maximum segment lengths
c
      if(cproc.ne.'1')then
        do j=1,nc
          f(j)=idate(j,3)
        enddo
        call mfs(f,nc,xmedseg,hhseg,hlseg,fsdseg,rmnseg,rmxseg)
        mnseg=int(rmnseg)
        mlseg=int(hlseg)
        medseg=int(xmedseg)
        mhseg=int(hhseg)
        mxseg=int(rmxseg)
        nsegmin=0
        nseghl=0
        nsegmed=0
        nseghh=0
        nsegmax=0
        do j=1,nc
          nyrs=idate(j,3)
          if(nyrs.eq.mnseg)nsegmin=nsegmin+1
          if(nyrs.ge.mlseg)nseghl=nseghl+1
          if(nyrs.ge.medseg)nsegmed=nsegmed+1
          if(nyrs.ge.mhseg)nseghh=nseghh+1
          if(nyrs.eq.mxseg)nsegmax=nsegmax+1
        enddo
        write(*,'(/1x,''CENSORED DATA SET'',t30,''   MIN   25%   MED   75%   MAX'')')
        write(*,'(1x,''Segment length statistics:'',t30,5i6)')mnseg,mlseg,medseg,mhseg,mxseg
        write(*,'(1x,''Series .ge. segment length:'',t30,5i6)')nsegmin,nseghl,nsegmed,nseghh,nsegmax
      endif
c
c *** get curve fit option
c
      write(*,'(/1x,''Curve fit options:''/
     *  1x,''Spline pct N filter (variable)  .le. -10''/
     *  1x,''Friedman variable span smoother'',t39,''-5''/
c    *  1x,''Regional curve detrending (RCS)'',t39,''-2''/
     *  1x,''Neg Expon/linear any slope'',t40,''1''/
     *  1x,''Neg Expon/linear neg slope'',t40,''2''/
     *  1x,''Linear regression any slope'',t40,''3''/
     *  1x,''Linear regression neg slope'',t40,''4''/
     *  1x,''Horizontal line through mean'',t40,''5''/
     *  1x,''Age-dependent spline smoothing'',t40,''8''/
     *  1x,''Median segment length spline'',t40,''9''/
     *  1x,''Spline filter in years (fixed)  .ge.  10'',t50,''==> '',$)')
      read(*,*,err=900)iopt
      if(iopt.eq.9)iopt=medseg
      if(iopt.eq.-5)then
        write(*,'(/1x,''Enter tweeter alpha 0-9'',t50,''==> '',$)')
        read(*,*)alpha
      endif
c
      switch=' '
      write(*,'(/1x,''Allow hierarchical curve switching? y/<n>'',t50,''==> '',$)')
      read(*,'(a)')switch
      if(switch.eq.'/')stop
      if(switch.eq.' ')switch='n'
      if(switch.eq.'Y')switch='y'
      if(switch.eq.'N')switch='n'
c
      itrun='n'
c      write(*,'(/1x,''Truncate curve fitting to same end year?''/
c     *  1x,''Curve is extrapolated to end. So y/<n>'',t50,''==> '',$)')
c      read(*,'(a)')itrun
c      if(itrun.eq.' ')itrun='n'
c      if(itrun.eq.'Y')itrun='y'
c      if(itrun.eq.'N')itrun='n'
c      if(itrun.eq.'y')then
c        write(*,'(1x,''Enter the fixed end year of fitting period'',t50,''==> '',$)')
c        read(*,*)itrunyr
c      endif
c
c      irev1=' '
c      write(*,'(/1x,''Review the curve fits? y/<n>'',t50,''==> '',$)')
c      read(*,'(a)')irev1
c      if(irev1.eq.' ')irev1='n'
c      if(irev1.eq.'Y')irev1='y'
c      if(irev1.eq.'N')irev1='n'
c      irev2=irev1
c
      irob1=' '
      write(*,'(/1x,''How to calculate the mean chronology:''/
     *           1x,''Non-robust arithmetic mean'',t39,''<1>''/
     *           1x,''Robust Tukey biweight mean'',t39,'' 2'',t50,''==> '',$)')
      read(*,'(a)')irob1
      if(irob1.eq.'/')stop
      if(irob1.eq.' ')irob=1
      if(irob1.eq.'1')irob=1
      if(irob1.eq.'2')irob=2
c
      reset=' '
      write(*,'(/1x,''Reset means of signal free series? <y>/n'',t50,''==> '',$)')
      read(*,'(a)')reset
      if(reset.eq.'/')stop
      if(reset.eq.' ')reset='y'
      if(reset.eq.'Y')reset='y'
      if(reset.eq.'N')reset='n'
c
      melvin=' '
      write(*,'(/1x,''H-P filter chronology differences? <y>/n'',t50,''==> '',$)')
      read(*,'(a)')melvin
      if(melvin.eq.'/')stop
      if(melvin.eq.' ')melvin='y'
      if(melvin.eq.'Y')melvin='y'
      if(melvin.eq.'N')melvin='n'
   10 itol=' '
      write(*,'(/1x,''Enter the signal free convergence test''/
     *           1x,''based on minimizing the mean absolute''/
     *           1x,''differences between the current and''/
     *           1x,''preceeding chronologies:  <1> 0.0005''/
     *           1x,''                           2  0.0001''/
     *           1x,''  *** roll your own ***    3  X.XXXX'',t50,''==> '',$)')
      read(*,'(a)')itol
      if(itol.eq.'/')stop
      if(itol.eq.' '.or.itol.eq.'1')then
        tol=0.0005d0
      elseif(itol.eq.'2')then
        tol=0.0001d0
      elseif(itol.eq.'3')then
        write(*,'(/1x,''Okay, enter YOUR convergence test'',t50,''==> '',$)')
        read(*,*,err=10)tol
      else
        goto 10
      endif
      write(*,'(/1x,''Okay, let us begin:''/)')
c
c      irevit=' '
c      write(*,'(/1x,''Review the signal free iterations? y/<n>'',t50,''==> '',$)')
c      read(*,'(a)')irevit
c      if(irevit.eq.' ')irevit='n'
c      if(irevit.eq.'Y')irevit='y'
c      if(irevit.eq.'N')irevit='n'
c
c **************************************************************************************************
c
c *** set the fixed truncation years if necessary
c
      do j=1,nc
        if(itrun.eq.'y')then
          ilyr=idate(j,2)
          if(itrunyr.lt.ilyr)then
            idate(j,3)=itrunyr
          else
            idate(j,3)=ilyr
          endif
        elseif(itrun.eq.'n')then
          idate(j,3)=idate(j,2)
        endif
      enddo
c
c *** do the signal free detrending
c
      call sigfree(x,y,f,ls,alpha,sfind,nc,crn1,crn2,nss,idate,ident,iopt,irob,
     *  switch,reset,melvin,tol,mxy,mxs,itrdb)
c
c *** open Tucson format chronology output file and save the standard and signal free chronologies
c
      call blank(fil1)
      fil1(1:40)=filnam(1:40)
      fil1(ll+1:ll+13)='_crns_tuc.txt'
      open(iu11,file=fil1,status='unknown')
      do j=1,3
        itrdb(j)(81:85)='  std'
        write(iu11,'(a)')itrdb(j)
      enddo
      call blank(id)
      id=filnam(1:6)
      call trrw(id,ifmny,maxn,crn1,nss,ifl,iu11,iuo,'iw',tit,'std',mxy,itrdb)
      do j=1,3
        itrdb(j)(81:85)='  ssf'
        write(iu11,'(a)')itrdb(j)
      enddo
      call trrw(id,ifmny,maxn,crn2,nss,ifl,iu11,iuo,'iw',tit,'ssf',mxy,itrdb)
      close(iu11)
c
c *** open column format chronology output file and save the standard and signal free chronologies
c
      fil1(ll+1:ll+13)='_crns_col.txt'
      open(iu11,file=fil1,status='unknown')
c      if(istab.eq.'n')then
      write(iu11,'(''  Year  SSize  StdCrn  SsfCrn  StdSig  SsfSig  StdSee  SsfSee'')')
c      else
c        write(iu11,'(''  Year  SSize  StdCrn  StdStb  SigFre  SFStab'')')
c      endif
      do i=1,maxn
        rns=float(nss(i))
        write(iu11,'(2i6,1x,10f8.3)')ifmny+i-1,nss(i),crn1(i,1),crn2(i,1),crn1(i,2),crn2(i,2),
     *    crn1(i,2)/sqrt(rns),crn2(i,2)/sqrt(rns)
      enddo
      close(iu11)
  900 pause
      stop
      end
*sbr blank: blanks character string
*
*  programmed dec 1986 by
*
*     richard  l  holmes
*     tree-ring laboratory
*     university of arizona
*     tucson, arizona  85721  usa
*
*  modified by r l holmes  9 dec 1986
*
      subroutine blank(l)
      character l*(*)
      ll=len(l)
      do 10 i=1,ll
   10 l(i:i)=' '
      return
      end
*sbr trread
*
      subroutine trread(x,y,nc,idate,ident,itrdb,iu10,mxy,mxs)
      real*4 x(mxy,mxs),y(mxy)
      integer idate(mxs,3)
      character ident(mxs)*8,id*8,itrdb(3)*85,tag*3,qq*2,tit*85
      ifmny=9999
      ilmxy=-9999
      ifl=0
      nc=0
    1 call blank(id)
      qq='mr'
      call trrw(id,iyr,n,y,nc,ifl,iu10,iuo,qq,tit,tag,mxy,itrdb)
      if(ifl.gt.0)then 
        nc=nc+1
        lyr=iyr+n-1
        idate(nc,1)=iyr
        idate(nc,2)=lyr
        idate(nc,3)=n
        ident(nc)=id
        ifmny=min0(iyr,ifmny)
        ilmxy=max0(iyr+n-1,ilmxy)
        do i=1,n
          x(i,nc)=y(i)
        enddo
        goto 1
      else
        idate(nc+1,1)=ifmny
        idate(nc+1,2)=ilmxy
        maxn=ilmxy-ifmny+1
        idate(nc+1,3)=maxn
      endif
      if(itrdb(1)(1:6).eq.'      ')then
        itrdb(1)(1:80)=tit(1:80)
        itrdb(1)(59:80)='1st itrdb line missing'
      endif
      if(itrdb(2)(1:6).eq.'      ')then
        itrdb(2)(1:49)=itrdb(1)(1:49)
        itrdb(2)(59:80)='2nd itrdb line missing'
      endif
      if(itrdb(3)(1:6).eq.'      ')then
        itrdb(3)(1:49)=itrdb(1)(1:49)
        itrdb(3)(59:80)='3rd itrdb line missing'
      endif
      return
      end
*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'
*      cru data,     write on file:    q = 'cw'
*
************************************************************
*
      subroutine trrw(id,iyr,n,y,nc,ifl,iu,iuo,q,tit,tag,mxy,itrdb)
      character id*8,ida*8,line*80,q*2,tit*85,tag*3
      character itrdb(3)*85
      integer nc(mxy),iy(10),inc(10),flag
      real*4 y(mxy)
      flag=-9999
ccc      tag=tit(1:3)
      call blank(tit)
      ntit=min0(len(tit),80)
      if(q .eq. 'ir')then
        k=1
        if(id(1:6) .eq. '******')k=-1
      elseif(q .eq. 'mw' .or. q .eq. 'iw' .or. q .eq. 'cr')then
        goto 5
      endif
    2 n=0
      isw=0
      ier=0
      ifl=max0(ifl,0)
      ifl=ifl+1
      if(ifl.eq.1)then
        do j=1,3
          call blank(itrdb(j))
        enddo
      endif
      goto 1
   91 ier=ier+1
      if(ier .gt. 25)stop'sbr trrw:  over 25 text lines'
      if(ier.eq.1)then
        if(iuo .gt. 0)write(iuo,'(/t10,''title of data:''/1x,a)')line
        write(*,'(/1x,''title of data:''/1x,a)')line
        tit(1:ntit)=line(1:ntit)
        itrdb(ier)(1:ntit)=line(1:ntit)
      else
        if(iuo .gt. 0)write(iuo,'(1x,a)')line
          if(ier.le.3)then
            itrdb(ier)(1:ntit)=line(1:ntit)
          endif
        write(*,'(1x,a)')line
      endif
    1 call blank(line)
      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,mxy)
        id=line(1:8)
        if(iuo .gt. 0)write(iuo,'(/t10,''title of data:  '',a/)')line
        return
      endif
      if(q .eq. 'ir')then
        if(line(7:10).eq.'    ')goto 1
        if(line(6:6).ne.'-'.and.line(5:5).ne.'-')then
          read(line,'(a6,t7,i4,10(i4,i3))',err=91,end=99)
     +      ida,iyr1,(iy(i),inc(i),i=1,10)
        elseif(line(6:6).eq.'-')then
          read(line,'(a5,t6,i5,10(i4,i3))',err=91,end=99)
     +      ida,iyr1,(iy(i),inc(i),i=1,10)
        elseif(line(5:5).eq.'-')then
          read(line,'(a4,t5,i6,10(i4,i3))',err=91,end=99)
     +      ida,iyr1,(iy(i),inc(i),i=1,10)
        endif
      ida(7:8)='  '
      elseif(q .eq. 'mr')then
        if(line(9:12).eq.'    ')goto 1
        if(line(8:8).ne.'-'.and.line(7:7).ne.'-')then
          read(line,'(a8,i4,10i6)',err=91,end=99)ida,iyr1,iy
        elseif(line(8:8).eq.'-')then
          read(line,'(a7,i5,10i6)',err=91,end=99)ida,iyr1,iy
        elseif(line(7:7).eq.'-')then
          read(line,'(a6,i6,10i6)',err=91,end=99)ida,iyr1,iy
        endif
      endif
ccc      write(*,'(1x,a,i8,10i6)')ida,iyr1,iy
      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
c
c *** write out tree-ring indices in tucson or cru format
c
    5 if(n.le.0)return
      ier=0
      j=0
      kyr=iyr
      if(q.eq.'iw'.or.q.eq.'cr')then
        if(q.eq.'cr')then
          ify=iyr
          ily=ify+n-1
          if(ify.le.-1000.or.ily.le.-1000)then
            write(iu,'(a6,2i5)')id(1:6),ify,ily
          else
            write(iu,'(a6,2i4)')id(1:6),ify,ily
          endif
        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   if(q.eq.'cr'.and.j.eq.n)return
        do i=1,10
          iy(i)=9990
          inc(i)=0
        enddo
        do 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
          iy(i)=nint(y(j)*1000.)
        enddo
   41   if(kyr.gt.-1000)then
          write(iu,'(a6,i4,10(i4,i3),2x,a)')
     +      id(1:6),kyr,(iy(k),inc(k),k=1,10),tag
        elseif(kyr.lt.-1000.and.kyr.gt.-10000)then
          write(iu,'(a5,i5,10(i4,i3),2x,a)')
     +      id(1:5),kyr,(iy(k),inc(k),k=1,10),tag
        else
          write(iu,'(a4,i6,10(i4,i3),2x,a)')
     +      id(1:5),kyr,(iy(k),inc(k),k=1,10),tag
        endif
        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.)
          if(iy(i) .eq. 999)iy(i)=1000
        endif      
      endif
   50 continue
   51 if(kyr.gt.-1000)then
        if(l.eq.10)then
          write(iu,'(a8,i4,10i6,t75,a)')id,kyr,(iy(k),k=1,l),tag
        else
          write(iu,'(a8,i4,10i6)')id,kyr,(iy(k),k=1,l)
        endif
      elseif(kyr.lt.-1000.and.kyr.gt.-10000)then
        if(l.eq.10)then
          write(iu,'(a7,i5,10i6,t75,a)')id,kyr,(iy(k),k=1,l),tag
        else
          write(iu,'(a7,i5,10i6)')id,kyr,(iy(k),k=1,l)
        endif
      else
        if(l.eq.10)then
          write(iu,'(a6,i6,10i6,t75,a)')id,kyr,(iy(k),k=1,l),tag
        else
          write(iu,'(a6,i6,10i6)')id,kyr,(iy(k),k=1,l)
        endif
      endif
      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,mxy)
      character line*80,frm*8,s*2,id*(*)
      real x(mxy)
      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
*sbr plot1
*
      subroutine plot1(nmax,wrk,sp,ny,iy,xlab,ylab,tit,rline,annot1,annot2)
c
c     write(winu9,'(/13x,''COLOR TABLE ______________________________________________'')')
c     write(winu9,'(13x,''White(0) Red(1)  Yellow(2)  Green(3)  Aquamarine(4) Pink(5) '')')
c     write(winu9,'(13x,''Wheat(6)  Grey(7)  Browm(8)  Blue(9)  BlueViolet(10)'')')
c     write(winu9,'(13x,''Cyan(11)  Turquoise(12)  Magenta(13)  Salmon(14)  Black(15)'')')
c     write(winu9,'(13x,''___________________________________________________________'',/)')
c
c
c     argument
      integer nmax,ny,iy
      dimension sp(nmax),wrk(nmax),rline(nmax)
      character tit*60,xlab*20,ylab*20,annot1*35,annot2*35
c     local
      integer p,red,green,blue
      real*4 x(10000),mny,mxy,mx,mnx	
      logical bg
c
      mny=100000.0;mxy=-100000.0
      red=255;green=255;blue=255
      bg=.false.
c
      do p=1,ny
        if(wrk(p).gt.mxy)mxy=wrk(p)
        if(sp(p).gt.mxy)mxy=sp(p)
        if(wrk(p).lt.mny)mny=wrk(p)
        if(sp(p).lt.mny)mny=sp(p)
        x(p)=iy+p-1
      enddo
      
      irange=mxy-mny
      mxy=mxy+irange/10
      mny=mny-irange/10
c      
      mnx=iy;mx=iy+ny-1
ccc      write(*,*)mxy,mny
ccc      write(*,*)mx,mnx
c      
      call plsdev('Mac1')
c
c set background color
c
      if(.not.bg)then
        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
      elseif(bg)then
        call plscol0(15,255,255,255)              !change default color map white>black
        call plscol0(0,0,0,0)                   !change default color map black>white
        call plscolbg(red,green,blue)            !change background color if other than white is desired
      endif  
c
c     initialize plplot
c
      call plinit()
c
c     set lable color
c
      call plcol(15)   !color box and lables
c     set curve width
      call plwid(2)
c
c     define the viewport one column by one rows
c
      call plssub(1,1)
ccc      call plenv(mnx,mx,mny,mxy,0,2)
      call plenv(mnx,mx,mny,mxy,0,0)
      if(rline(1).ne.-999.9)then
        call plline(ny,x,rline)            !plot rline
      endif
      call plfont(1)
      call pllab(xlab,ylab,tit)
c     set line style
      call pllsty(1)
c     set line color
      call plcol(15)  !color lines
      call plline(ny,x,wrk)            !plot wrk
ccc      call plmtex('b',-25.5,0.34,0.0,' Original')
      call plmtex('b',-29.0,0.04,0.0,annot1)
c     set curve width
      call plwid(2)
c     set curve color
      call plcol(1)  !color lines
      call plline(ny,x,sp)           !plot spline curve
ccc      call plmtex('b',-29.5,0.52,0.0,' Signal Free')
      call plmtex('b',-27.0,0.04,0.0,annot2)
c     set curve width
ccc      call plwid(2)
c     set lable color
ccc      call plcol(15)   !color box and lables
ccc      call plline(ny,x,cl)           !plot mean curve
c     
      call plend()
c     
      return
      end
*sbr spaghetti
*
      subroutine spaghetti(xy,ny,ns,ify,xlab,ylab,tit,annot1,annot2,mxy,mxs)
c
c     write(winu9,'(/13x,''COLOR TABLE ______________________________________________'')')
c     write(winu9,'(13x,''White(0) Red(1)  Yellow(2)  Green(3)  Aquamarine(4) Pink(5) '')')
c     write(winu9,'(13x,''Wheat(6)  Grey(7)  Browm(8)  Blue(9)  BlueViolet(10)'')')
c     write(winu9,'(13x,''Cyan(11)  Turquoise(12)  Magenta(13)  Salmon(14)  Black(15)'')')
c     write(winu9,'(13x,''___________________________________________________________'',/)')
c
      parameter(mxlen=10000)
c     argument
      real*4 xy(mxy,mxs)
      integer ny,ns,ify,mxy,mxs
      character tit*60,xlab*20,ylab*20,annot1*35,annot2*35
c     local
      integer red,green,blue
      real*4 x(mxlen),wrk(mxlen),miny,maxy,maxx,minx,range	
      logical bg
c
      mny=100000.0;mxy=-100000.0
      red=255;green=255;blue=255
      bg=.false.
c
      minx=ify
      maxx=ify+ny-1
      maxy=-100000.0
      miny=100000.0
      do i=1,ny
        do j=1,ns
          datum=xy(i,j)
          if(datum.ne.-9.99)then
            if(datum.gt.maxy)maxy=datum
            if(datum.lt.miny)miny=datum
          endif
        enddo
        x(i)=ify+i-1
      enddo
      range=maxy-miny
      maxy=maxy+range/10.
      miny=miny-range/10.
c      
      call plsdev('Mac1')
c
c set background color
c
      if(.not.bg)then
        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
      elseif(bg)then
        call plscol0(15,255,255,255)              !change default color map white>black
        call plscol0(0,0,0,0)                   !change default color map black>white
        call plscolbg(red,green,blue)            !change background color if other than white is desired
      endif  
c
c *** initialize plplot
c
      call plinit()
c
c *** set lable color
c
      call plcol(15)   !color box and lables
c
c *** set curve width
c
      call plwid(2)
c
c *** define the viewport one column by one rows
c
      call plssub(1,1)
      call plenv(minx,maxx,miny,maxy,0,0)
      call plfont(1)
      call pllab(xlab,ylab,tit)
c
c *** set line style
c
      call pllsty(1)
c
c *** set line color
c
      call plcol(15)  !color lines
c
c *** do plots
c
      do j=1,ns
        do i=1,ny
          wrk(i)=xy(i,j)
        enddo
        call plwid(1)
        call plcol(15)
        if(j.eq.1)then
          call plcol(1)
          call plwid(2)
        elseif(j.eq.ns)then
          call plcol(9)
          call plwid(2)
        endif
        call plline(ny,x,wrk)
      enddo
      call plwid(2)
      call plcol(1)
      call plmtex('b',-29.0,0.05,0.0,annot1)
      call plcol(9)
      call plmtex('b',-27.0,0.05,0.0,annot2)
c
c *** end plotting
c
      call plend()     
      return
      end
*sbr sigfree
*
      subroutine sigfree(x,y,f,ls,alpha,sfind,nc,crn1,crn2,nss,idate,ident,iopt,irob,
     *  switch,reset,melvin,tol,mxy,mxs,itrdb)
      parameter(mxlen=10000,mxser=1000,maxit=15,iu29=29,iu30=30,iu31=31,iu32=32)
      real*4 x(mxy,mxs),y(mxy),f(mxy),sfind(mxy,mxs),crn1(mxy,2),crn2(mxy,2),savit(mxlen,maxit)
      real*4 crnsav(mxlen),xbar(mxser),segmn(mxlen),segsd(mxlen),agemn(mxlen),agesd(mxlen)
      real*4 test1(maxit),test2(maxit),test3(maxit),crns(mxlen,maxit+1),wrk1(mxlen),wrk2(mxlen)
      real*8 test,tol,testd,rough
      integer idate(mxs,3),nss(mxy),iopt1(mxser,maxit+2),ls(mxy)
      character ident(mxs)*8,tit*60,xlab*20,ylab*20,annot1*35,annot2*35,melvin,line*410,reset,switch
      character star1,star2,star3,itrdb(3)*85,id*8
c
c *** get chronologiy length info
c
      ifmny=idate(nc+1,1)
      ilmxy=idate(nc+1,2)
      maxn=ilmxy-ifmny+1
c
c *** calculate the average segment length over all years
c
      do j=1,nc
        ify=idate(j,1)
        ily=idate(j,2)
        n=ily-ify+1
        ish=ify-ifmny
        do i=1,n
          ii=i+ish
          sfind(ii,j)=float(n)
        enddo
      enddo
      call robsum(sfind,maxn,nc,segmn,segsd,nss,1,mxy,mxs)
c
c *** calculate the average tree age over all years
c
      do j=1,nc
        ify=idate(j,1)
        ily=idate(j,2)
        n=ily-ify+1
        ish=ify-ifmny
        do i=1,n
          ii=i+ish
          sfind(ii,j)=float(i)
        enddo
      enddo
      call robsum(sfind,maxn,nc,agemn,agesd,nss,1,mxy,mxs)
c
c *** plot out mean segment length and mean tree age
c
      call blank(tit)
      tit(1:32)='Mean Segment Length And Mean Age'
      call blank(xlab)
      xlab(1:5)='Years'
      call blank(ylab)
      ylab(1:5)='Years'
      call blank(annot1)
      annot1(1:8)='Mean Seg'
      call blank(annot2)
      annot2(1:8)='Mean Age'
      do i=1,maxn
        f(i)=0.0
      enddo
      call plot1(mxy,segmn,agemn,maxn,ifmny,xlab,ylab,tit,f,annot1,annot2)
c
c *** do the initial detrending based on iopt curve fitting option
c
      open(iu29,file='Chronologies_Iter=0-Final.txt',status='unknown')
      close(iu29,status='delete')
      open(iu29,file='Chronologies_Iter=0-Final.txt',status='unknown')

      open(iu31,file='OriginalCurves.txt',status='unknown')
      close(iu31,status='delete')
      open(iu31,file='OriginalCurves.txt',status='unknown')
      open(iu32,file='OriginalIndices.txt',status='unknown')
      close(iu32,status='delete')
      open(iu32,file='OriginalIndices.txt',status='unknown')
      do j=1,nc
        ify=idate(j,1)
        ily=idate(j,2)
        ily1=ily
        n=ily-ify+1
        n1=n
        do i=1,n
          y(i)=x(i,j)
        enddo
        call mnsd(y,n,xm,sd)
        xbar(j)=xm
        iopt1(j,1)=iopt
        iopt1(j,2)=iopt
        call curvefits(y,f,ify,n,n1,tit,iopt1(j,2),alpha,mxy,mxs)
        do i=1,n
          wrk2(i)=y(i)/f(i)
          sfind(i,j)=y(i)/f(i)
        enddo

        call trrw(ident(j),ify,n,f,nss,ifl,iu31,iuo,'mw',tit,'crv',mxy,itrdb)
        call trrw(ident(j),ify,n,wrk2,nss,ifl,iu32,iuo,'mw',tit,'ind',mxy,itrdb)

      enddo
      close(iu31)
      close(iu32)
c
c *** align the data by calendar year now
c
      do j=1,nc
        ify=idate(j,1)
        ily=idate(j,2)
        n=ily-ify+1
        ish=ify-ifmny
        if(ish.gt.0)then
          do i=1,n
            ii=n-i+1
            x(ii+ish,j)=x(ii,j)
            sfind(ii+ish,j)=sfind(ii,j)
          enddo
          do i=1,ish
            x(i,j)=-9.99
            sfind(i,j)=-9.99
          enddo
        endif
      enddo
c
c *** calculate the initial mean chronology
c
      call robsum(sfind,maxn,nc,crn1(1,1),crn1(1,2),nss,irob,mxy,mxs)

      nss1=0
      do i=1,maxn
        crn2(i,1)=crn1(i,1)
        crnsav(i)=crn1(i,1)
        crns(i,1)=crnsav(i)
      enddo
      do i=1,maxn
        if(nss(i).eq.1)then
          nss1=nss1+1
        else
          goto 5
        endif
      enddo
    5 if(nss1.gt.0)then
        write(*,'(/1x,''nss1 skip first number of years: '',i3)')nss1
        do j=1,nc
          ify=idate(j,1)
          ily=idate(j,2)
          n=ily-ify+1
          if(ify.eq.ifmny)then
            ncset=j
            nreset=n-nss1
            write(*,'(/1x,i4,2x,a8,2x,5i6)')ncset,ident(j),ify,ily,n,nss1,nreset
            do i=1,nreset
              ii=i+nss1
              y(i)=x(ii,j)
            enddo
            call mnsd(y,nreset,xm1,sd)
            write(*,'(/1x,''full and truncated means: '',2f8.4)')xbar(j),xm1
            xbar(j)=xm1
          endif
        enddo
      endif

      id='iter00'
      call trrw(id,ifmny,maxn,crns(1,1),nss,ifl,iu29,iuo,'iw',tit,'crn',mxy,itrdb)

c
c *** do the signal free iterations
c
      iter=0
      small=9999.0
      smooth=9999.0
      big=9999.0
      write(*,'(38x,''Mean ABS      Mean dABS      Maxx ABS'')')
   10 continue
      iter=iter+1
      open(iu30,file='SigFreeRawdat.txt',status='unknown')
      close(iu30,status='delete')
      open(iu30,file='SigFreeRawdat.txt',status='unknown')
      open(iu31,file='SigFreeCurves.txt',status='unknown')
      close(iu31,status='delete')
      open(iu31,file='SigFreeCurves.txt',status='unknown')
      open(iu32,file='SigFreeIndices.txt',status='unknown')
      close(iu32,status='delete')
      open(iu32,file='SigFreeIndices.txt',status='unknown')

      do j=1,nc
        ify=idate(j,1)
        ily=idate(j,2)
        ily1=idate(j,3)
        n=ily-ify+1
        n1=ily1-ify+1
        iskp=ify-ifmny
        do i=1,n
          ii=i+iskp
          if(nss(ii).gt.1)then
            y(i)=x(ii,j)/crn2(ii,1)
          else
            y(i)=x(ii,j)
          endif
        enddo
        
        if(reset.eq.'y')then
          if(nss1.gt.0.and.j.eq.ncset)then
            call mnsd(y(nss1+1),nreset,xm1,sd)
          else
            call mnsd(y,n,xm1,sd)
          endif
          do i=1,n
            y(i)=y(i)/xm1*xbar(j)
          enddo
        endif

ccc        if(reset.eq.'y')then
ccc          call mnsd(y,n,xm1,sd)
ccc          do i=1,n
ccc            y(i)=y(i)/xm1*xbar(j)
ccc          enddo
ccc        endif
        
        iopt1(j,iter+2)=iopt
        call curvefits(y,f,ify,n,n1,id,iopt1(j,iter+2),alpha,mxy,mxs) 
        if(switch.eq.'n'.and.iter.ge.5)then
          if((iopt1(j,iter+2).ne.iopt1(j,iter+1)).or.(iopt1(j,iter+1).eq.-66))then
            iopt1(j,iter+2)=-66
            call curvefits(y,f,ify,n,n1,id,iopt1(j,iter+2),alpha,mxy,mxs)      
          endif
        endif
        do i=1,n
          ii=i+iskp
          wrk1(i)=y(i)
          y(i)=x(ii,j)
          sfind(ii,j)=y(i)/f(i)
          wrk2(i)=y(i)/f(i)
        enddo        

        call trrw(ident(j),ify,n,wrk1,nss,ifl,iu30,iuo,'mw',tit,'sfr',mxy,itrdb)
        call trrw(ident(j),ify,n,f,nss,ifl,iu31,iuo,'mw',tit,'sfc',mxy,itrdb)
        call trrw(ident(j),ify,n,wrk2,nss,ifl,iu32,iuo,'mw',tit,'sfi',mxy,itrdb)

      enddo
      close(iu30)
      close(iu31)
      close(iu32)
c
c *** now calculate the signal free mean chronology
c
      call robsum(sfind,maxn,nc,crn2(1,1),crn2(1,2),nss,irob,mxy,mxs)
      
      do i=1,maxn
        crns(i,iter+1)=crn2(i,1)
      enddo

      write(id(5:6),'(i2)')iter
      if(iter.lt.10)id(5:5)='0'
      call trrw(id,ifmny,maxn,crns(1,iter+1),nss,ifl,iu29,iuo,'iw',tit,'crn',mxy,itrdb)

      do i=1,maxn
        y(i)=crnsav(i)-crn2(i,1)
        savit(i,iter)=y(i)
      enddo
c
c *** high-pass filter the differences according to Tom's average segment length spline smoothing
c
      if(melvin.eq.'y')then
        do i=1,maxn-2
          ls(i)=segmn(i+1)
        enddo
        call spline2(maxn,y,f,ls)
        do i=1,maxn
          y(i)=y(i)-f(i)
        enddo
      endif
c
c *** compute the convergence statistics - use Tom's average absolute difference criterion
c
      do i=1,maxn
        f(i)=abs(y(i))
      enddo
      test=0.0d0
      testd=0.0d0
      maxn1=0
      do i=1,maxn
        if(nss(i).gt.1)then
          maxn1=maxn1+1
          f(maxn1)=f(i)
          test=test+dabs(y(i))
          if(i.gt.1)then
            testd=testd+dabs(y(i)-y(i-1))
          endif
        endif
      enddo
      test=test/dfloat(maxn1)
      if(test.lt.small)then
        small=test
        star1='*'
      else
        star1=' '
      endif
      testd=testd/dfloat(maxn1-1)
      if(testd.lt.smooth)then
        smooth=testd
        star2='*'
      else
        star2=' '
      endif
      call qcksrt(maxn1,f)
      testb=f(maxn1)
      if(testb.lt.big)then
        big=testb
        star3='*'
      else
        star3=' '
      endif
      test1(iter)=test
      test2(iter)=testd
      test3(iter)=testb
      write(*,'(1x,''Iter '',i3,''  Convergence criteria: '',f13.8,a,f14.8,a,f13.8,a)')
     *  iter,test,star1,testd,star2,testb,star3
c
c *** test for convergence now using the user-supplied mean absolute difference threshold
c *** but use the mean first-difference of the absolute differences as an additional test
c
      if(test.le.tol)then
        goto 20
      elseif(iter.eq.maxit)then
        goto 20
c      elseif(iter.ge.10.and.(test1(iter).gt.test1(iter-1)).and.(test2(iter).gt.test2(iter-1)))then
c        do i=1,maxn
c          crn2(i,1)=crnsav(i)
c        enddo
c        iter=iter-1
c        goto 20
      elseif(test.gt.tol)then
        do i=1,maxn
          crnsav(i)=crn2(i,1)
        enddo
        goto 10
      endif
c
c *** plot out the original and signal free chronologies
c
   20 write(*,'(/1x,''Final iteration for signal free: '',i2)')iter
      call blank(tit)
      if(irob.eq.1)then
        tit(1:41)='Original & Signal Free Chronology, Iter ='
        write(tit(42:44),'(i3)')iter
      elseif(irob.eq.2)then
        tit(1:48)='Robust Original & Signal Free Chronology, Iter ='
        write(tit(49:51),'(i3)')iter
      endif
      call blank(ylab)
      ylab(1:7)='Indices'
      call blank(annot1)
      annot1(1:8)='Original'
      call blank(annot2)
      annot2(1:11)='Signal Free'
      do i=1,maxn
        f(i)=1.0
      enddo
      call plot1(mxy,crn1(1,1),crn2(1,1),maxn,ifmny,xlab,ylab,tit,f,annot1,annot2)
c
c *** plot out the mean abs and mean dabs convergence statistics
c
      call blank(tit)
      tit(1:40)='Iteration Convergence Statistics, Iter ='
      write(tit(41:43),'(i3)')iter
      call blank(xlab)
      xlab(1:10)='Iterations'
      call blank(ylab)
      ylab(1:17)='Differences'
      call blank(annot1)
      annot1(1:14)='Mean  ABS Diff'
      call blank(annot2)
      annot2(1:14)='Mean dABS Diff'
      do i=1,maxn
        f(i)=0.0
      enddo
      call plot1(maxit,test1,test2,iter,1,xlab,ylab,tit,f,annot1,annot2)
c
c *** plot out the successive signal free iteration differences
c
      call blank(tit)
      tit(1:40)='Signal Free Iteration Differences'
      call blank(xlab)
      xlab(1:5)='Years'
      call blank(ylab)
      ylab(1:17)='Differences'
      call blank(annot1)
      annot1(1:15)='First Iteration'
      call blank(annot2)
      annot2(1:15)='Final Iteration'
      do j=2,iter
        jj=iter-j+1
        rj=j-1
        do i=1,maxn
          savit(i,jj)=savit(i,jj)-(rj*0.10)
        enddo
      enddo
      call spaghetti(savit,maxn,iter,ifmny,xlab,ylab,tit,annot1,annot2,mxlen,maxit) 
c
c *** print out the curve fit options as a function of signal free iteration
c
      write(*,'(/1x,''SIGNAL FREE DETRENDING OPTIONS BY ITERATION:'')')
      write(*,'(3x,''NC  IDENT'',4x,''   0    0'',100i5)')(j,j=1,iter)
      call blank(line)
      do i=1,iter+2
        ii=(i-1)*5+1
        iii=ii+4
        line(ii:iii)=' ---'
      enddo
      write(*,'(3x,''--  -----'',4x,a)')line
      do j=1,nc
        write(*,'(i5,2x,a,102i5)')j,ident(j),(iopt1(j,k),k=1,iter+2)
      enddo
      close(iu29)
      return
      end
*sbr curvefits
*
      subroutine curvefits(y,f,iyr,n,n1,id,iopt,alpha,mxy,mxs)
      parameter(mxlen=10000,epsi=0.000001)
      real*4 y(mxy),f(mxy),eq(7)
      real*4 wrk(mxlen),w(mxlen),sc(mxlen,7),tmp(mxlen)
      character id*8
c
c *** negative exponential curve
c
      iopt1=iopt
      if(iopt.eq.1.or.iopt.eq.2)then
        do j=1,7
          eq(j)=0.0
        enddo
        call curve(n,n1,y,f,eq)
        if(eq(2).le.epsi.or.eq(5).lt.epsi)then
ccc        if(eq(2).le.0.or.eq(5).lt.0)then
          if(iopt.eq.1)iopt=3
          if(iopt.eq.2)iopt=4
          do j=1,7
            eq(j)=0.0
          enddo
          call trend(n,n1,y,f,eq(3),eq(5))
        endif
        if(eq(3).gt.epsi.and.iopt1.eq.2)then
ccc        if(eq(3).gt.0.and.iopt1.eq.2)then
          iopt=5
          call xmean(y,f,n,n1)
        endif
c
c *** linear regression curve
c
      elseif(iopt.eq.3.or.iopt.eq.4)then
        do j=1,7
          eq(j)=0.0
        enddo
        call trend(n,n1,y,f,eq(3),eq(5)) 
        if(iopt.eq.4.and.eq(3).gt.0)then
          iopt=5
          call xmean(y,f,n,n1)
        endif
c
c *** horizontal line through the mean
c
      elseif(iopt.eq.5)then
        call xmean(y,f,n,n1)
c
c *** Friedman variable span smoother
c
      elseif(iopt.eq.-5)then
        do i=1,n
          wrk(i)=float(i)
          w(i)=1.0
        enddo
        span=0.0
        if(n1.lt.40)span=0.3
        iper=1
        call supsmu(n1,wrk,y,w,iper,span,alpha,f,sc)
        if(n.gt.n1)then
          do i=1,n1
            tmp(i)=y(i)
          enddo
          iext=n-n1
          do i=1,iext
            ii=i+n1
ccc            f(ii)=f(n1)
            tmp(ii)=f(n1)
          enddo
          call supsmu(n,wrk,tmp,w,iper,span,alpha,f,sc)
        endif
c
c *** cubic smoothing spline smoothing
c
      elseif(abs(iopt).ge.8)then
        if(iopt.eq.8)then
          ls=iopt
          call spline1(n1,y,f,ls)
          
c          if(n.eq.n1)then
c            minspl=n
c            do i=2,n
c              ii=n-i+1
c              iii=n-i
c              dif=f(ii)-f(iii)
c              if(dif.le.0)then
c                minspl=n-ii
c                goto 10
c              endif
c            enddo
c   10       if(minspl.lt.n)then
c              do i=minspl,n
c                f(i)=f(minspl-1)
c              enddo
c            endif
c          endif
          
        elseif(iopt.ge.10)then
          ls=iopt
          call spline1(n1,y,f,ls)
        elseif(iopt.le.-10)then
          rn1=n1
          ropt=float(abs(iopt))/100.
          ls=ropt*rn1
          call spline1(n1,y,f,ls)
        endif
        if(n.gt.n1)then
          do i=1,n1
            tmp(i)=y(i)
          enddo
          iext=n-n1
          do i=1,iext
            ii=i+n1
ccc            f(ii)=f(n1)
            tmp(ii)=f(n1)
          enddo
          call spline1(n,tmp,f,ls)
        endif
      endif
      
      ier=0
      do i=1,n
        if(f(i).lt.epsi)then
          ier=1
          goto 90
        endif
      enddo
   90 if(ier.eq.1)then
        iopt=5
        call xmean(y,f,n,n1)
        goto 100
      endif
      do i=1,n
        if(f(i).lt.0.02)then
          f(i)=0.02
        endif
      enddo
      
c      ier=0
c      do i=1,n
c        if(f(i).le.0.02)then
c          ier=ier+1
c          goto 90
c        endif
c      enddo
c   90 if(ier.gt.0)then
c        if(ier.eq.1)iopt=-66
c        if(ier.eq.2)iopt=-50
c        if(ier.eq.3)iopt=-33
c        if(ier.eq.4)iopt=5
c        rn1=n1
c        ropt=float(abs(iopt))/100.
c        ls=ropt*rn1
c        if(iopt.lt.4)then
c          call spline1(n1,y,f,ls)
c        else
c          call xmean(y,f,n,n1)
c        endif
c
c        do i=1,n
c          if(f(i).le.0.02)then
c            ier=ier+1
c            goto 90
c          endif
c        enddo

ccc        call xmean(y,f,n,n1)
ccc        call plot1(mxy,y,f,n,iyr,id)
c        do i=1,n1
c          wrk(i)=float(i)
c          w(i)=1.0
c        enddo
c        span=0.0
c        if(n1.lt.40)span=0.3
c        iper=1
c        call supsmu(n1,wrk,y,w,iper,span,5,f,sc) 
c        iopt=-5
c        if(n.gt.n1)then
c          iext=n-n1
c          do i=1,iext
c            ii=i+n1
c            f(ii)=f(n1)
c          enddo
c        endif
ccc        call plot1(mxy,y,f,n,iyr,id)
c      endif
  100 return
      end
*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
*
**************************************************
*
      subroutine curve(n,n1,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,rn1,dlt,ab,arg,zz,dinc
      data dinc/.0000005d0/
      nit=0
      if(n1 .lt. 20)then
      write(*,'(''>> sbr curve:  only'',i3,
     +'' values, no exponential curve fit'')')n
      goto 99
      endif
      rn1=dble(n1)
      yf=0.d0
      do 5 i=1,10
    5 yf=yf+dble(z(i))
      yf=yf/10.d0
      yl=0.d0
      do 8 i=n1-9,n1
    8 yl=yl+dble(z(i))
      yl=yl/10.d0
      bh(2)=(dlog(yf)-dlog(yl))/(rn1-10.d0)
      if(bh(2) .le. 0.d0)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,n1
      arg=ab*dfloat(i)
      if(arg .gt. 112.8d0)goto 99
   19 f(i)=dexp(-bh(k)*dfloat(i))
      a=dble(n1)
      b=0.d0
      d=0.d0
      e=0.d0
      g=0.d0
      do 21 i=1,n1
      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,n1
      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
      do i=1,n							!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
*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
*
********************************************************************
*
      subroutine trend(n,n1,y,f,sl,yi)			!Hakan
      implicit real*8 (a-h,o-z)
      real*4 y(n),f(n),sl,yi
      an1=dfloat(n1)
      t=0.d0
      w=0.d0
      x=0.d0
      z=0.d0
      do i=1,n1
        ai=dfloat(i)
        z=z+ai
        t=t+ai*ai
        x=x+dble(y(i))*ai
        w=w+dble(y(i))
      enddo
      ymean=w/an1
      dsl=(x-w*z/an1)/(t-z*z/an1)
      dyi=ymean-dsl*z/an1
      do i=1,n						!Hakan
        ai=dfloat(i)
        f(i)=sngl(dsl*ai+dyi)
      enddo
      sl=sngl(dsl)
      yi=sngl(dyi)
      return
      end
*sbr xmean
*
      subroutine xmean(x,f,n,n1)
      real*4 x(n),f(n)
      if(n1.le.0)then
        write(*,'(/1x,''error in xmean: n=0'')')
        pause
        stop
      endif
      rn1=n1
      xm=0.0
      do i=1,n1
        xm=xm+x(i)
      enddo
      xm=xm/rn1
      do i=1,n
        f(i)=xm
      enddo
      return
      end
      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
*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 filter 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: stiffness (50% frequency response) of spline
* a:  working array of dimension (n,4)
*
********************************************************************
*
      subroutine spline(n,y,f,pp,ls)
      parameter(mxy=10000)
      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=dble(ls)
      arg=(2.d0*pi)/v
      p=(6.d0*(dcos(arg)-1.d0)**2.d0)/(dcos(arg)+2.d0)
      pp=sngl(p)
      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
*sbr spline1: 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 filter 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: stiffness (50% frequency response) of spline
* a:  working array of dimension (n,4)
*
********************************************************************
*
      subroutine spline1(n,y,f,ls)
      parameter(mxy=10000)
      implicit real*8 (a-h,o-z)
      dimension a(mxy,4),c1(4),c2(3),p(mxy)
      real*4 y(n),f(n)
      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
      if(ls.eq.9)then
        ls1=20
      else
        ls1=ls
      endif
      
ccc      rn23=dfloat(n)*2.d0/3.d0
      rn23=dfloat(n)
      do i=1,nm2
        if(ls.eq.9)then
          v=dfloat(ls1)+i-1
          if(v.gt.rn23)v=rn23
        else
          v=dfloat(ls1)
        endif
        arg=(2.d0*pi)/v
        p(i)=(6.d0*(dcos(arg)-1.d0)**2.d0)/(dcos(arg)+2.d0)
      enddo
      
      do 1 i=1,nm2
      do 1 j=1,3
      a(i,j)=c1(j)+p(i)*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
*sbr spline2: 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 filter 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: stiffness (50% frequency response) of spline
* a:  working array of dimension (n,4)
*
********************************************************************
*
      subroutine spline2(n,y,f,ls)
      parameter(mxy=10000)
      implicit real*8 (a-h,o-z)
      integer ls(mxy)
      dimension a(mxy,4),c1(4),c2(3),p(mxy)
      real*4 y(n),f(n)
      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
      do i=1,nm2
        v=dfloat(ls(i))
        arg=(2.d0*pi)/v
        p(i)=(6.d0*(dcos(arg)-1.d0)**2.d0)/(dcos(arg)+2.d0)
      enddo
      
      do 1 i=1,nm2
      do 1 j=1,3
      a(i,j)=c1(j)+p(i)*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
*sbr robsum
*
      subroutine robsum(x,n,nc,crnmn,crnsd,nss,iopt,mxy,mxs)
      parameter(maxnum=1000)
      real*4 x(mxy,mxs),crnmn(mxy),crnsd(mxy),wrk(maxnum),datum
      integer nss(mxy),iopt
      if(nc.le.3)iopt=2
      do i=1,n
        n1=0
        do j=1,nc
          datum=0.0
          datum=x(i,j)
          if(datum.ne.-9.99)then
            n1=n1+1
            wrk(n1)=datum
          endif
        enddo
        if(n1.eq.1)then
          xm=wrk(1)
          sd=0.0
        elseif(n1.gt.1.and.n1.lt.6)then
          if(iopt.eq.1)then
            call mnsd(wrk,n1,xm,sd)
          elseif(iopt.eq.2)then
            call mdian1 (wrk,n1,xmed)
            call xmad(wrk,n1,xmed,xmadm)
            xm=xmed
            sd=3.0*xmadm/2.0
          endif
        elseif(n1.ge.6)then
          if(iopt.eq.1)then
            call mnsd(wrk,n1,xm,sd)
          elseif(iopt.eq.2)then
            call xbiwt(wrk,n1,xlbiwt,xsbiwt,xlbiwt1,xsbiwt1)
            xm=xlbiwt
            sd=xsbiwt
          endif
        endif
        crnmn(i)=xm
        crnsd(i)=sd
        nss(i)=n1
      enddo
      return
      end
*sbr mnsd
*
      subroutine mnsd(x,n,xm,sd)
      real*4 x(n),xm,sd,sq,sm,rn
      xm=0.0
      sd=0.0
c
c *** compute the mean and standard deviation of the data
c
      sq=0.0
      sm=0.0
      do 10 i=1,n
        sm=sm+x(i)
   10 continue
      rn=n
      xm=sm/rn
      do 15 i=1,n
        sq=sq+(x(i)-xm)**2
   15 continue
      sd=sqrt(sq/(rn-1.0))
      return
      end
c-----------------------------------------------------------------------------
      SUBROUTINE XBIWT (XDATA,N,XLBIWT,XSBIWT,XLBIWT1,XSBIWT1)
c-----------------------------------------------------------------------------
c
c     The subroutine XBIWT provides an estimator of the location and
c     scale of the data set XDATA.  The scale uses the Biweight function
c     in the general formula of "A-estimators." This formula is given
c     on page of 416 in UREDA (formula 4). The BIWEIGHT scale estimate
c     is returned as the value XSBIWT. The BIWEIGHT function is given
c     by:
c
c          u((1-u*u)**2)     abs(u) <= 1
c          f(u) =
c                                  0                 abs(u) >  1
c
c     where u is defined by
c
c          u = (XDATA(I) - M) / c*MAD  .
c
c     M, MAD, and c are the median, the median absolute deviation from
c     the median, and the tuning constant respectively. The tuning
c     constant is a parameter which is chosen depending on the sample
c     size and the specific function being used for the scale estimate.
c     (See page 417 in UREDA).  Here we take c = 9.0.
c
c     The biweght location is found using the formula:
c
c          T = M + (sums)
c
c          where M is the sample median and sums are
c          as given on page 421 in UREDA
c
c          the tuning constant c is set to 6.0 for calculation
c          of the location as reccommended by Tukey ()
c
c     NOTE that the biweight is meant to be an iterated estimator, but one
c     commonly only takes the first step of the iteration.  Here we report
c     both the one-step estimators (XLBIWT1, XSBIWT1) and the preferred
c     fully iterated versions (XLBIWT, XSBIWT).
c
c****************************************************************************
      parameter(iii=10000)
      implicit real*4 (a-h,o-z)
      dimension xdata(n),u1(iii),u2(iii),xlb(11),xsb(11)
      data zero,d6,d1,d5,d9/0.0,6.0,1.0,5.0,9.0/
 
c---  sort the data and find the median
 
      CALL MDIAN1(XDATA,N,XM)
 
c---  call xmad to find the median absolute deviation
 
      CALL XMAD(XDATA,N,XM,XMADM)
 
c---  must choose value of the tuning constant "c"
c     here c = 6.0 for the location estimator and
c     9.0 for the scale estimator
 
      c1 = d6
      c2 = d9
      
c      c1 = 7.5			!only for comparing with Tom's signal free
c      c2 = 7.5			!only for comparing with Tom's signal free
 
      if (xmadm.le..0001) then
        xlbiwt=xm
        xlbiwt1=xm
        xsbiwt=xmadm
        xsbiwt1=xmadm
        goto 20
      endif
 
      do i = 1,n
        u1(i) = (xdata(i) - xm)/(c1*xmadm)
        u2(i) = (xdata(i) - xm)/(c2*xmadm)
      enddo
      s1 = zero
      s2 = zero
      s3 = zero
      s4 = zero
      do i = 1,n
        if (abs(u2(i)) .lt. d1) then
          s1 = s1+(((xdata(i)-xm)**2)*(d1-(u2(i)*u2(i)))**4)
          s2 = s2+((d1-u2(i)*u2(i))*(d1-(d5*u2(i)*u2(i))))
        endif
        if (abs(u1(i)) .lt. d1) then
          s3 = s3+(xdata(i)-xm)*(d1-u1(i)*u1(i))**2
          s4 = s4+(d1-u1(i)*u1(i))**2
        endif
      enddo 
c---  here are the one-step estimators
      xlbiwt1 = xm+s3/s4
      xsbiwt1 = float(n)/(float(n-1))**0.5*s1**0.5/abs(s2)
c---  now obtain the fully-iterated versions
c---  solve for new estimates of u1 and u2
      xlb(1) = xlbiwt1
      xsb(1) = xsbiwt1
      do j = 2,11   !assuming 10 iterations is sufficient
        xmm = xlb(j-1)
        do i = 1,n
          u1(i) = (xdata(i) - xmm)/(c1*xmadm)
          u2(i) = (xdata(i) - xmm)/(c2*xmadm)
        enddo
        s1 = zero
        s2 = zero
        s3 = zero
        s4 = zero
        do i = 1,n
          if (abs(u2(i)) .lt. d1) then
            s1 = s1+(((xdata(i)-xmm)**2)*(d1-(u2(i)*u2(i)))**4)
            s2 = s2+((d1-u2(i)*u2(i))*(d1-(d5*u2(i)*u2(i))))
          endif
          if (abs(u1(i)) .lt. d1) then
            s3 = s3+(xdata(i)-xmm)*(d1-u1(i)*u1(i))**2
            s4 = s4+(d1-u1(i)*u1(i))**2
          endif
        enddo
        xlb(j) = xlb(j-1)+s3/s4
        xsb(j) = float(n)/(float(n-1))**0.5*s1**0.5/abs(s2)
      enddo
      xlbiwt = xlb(11)
      xsbiwt = xsb(11)
   20 continue
      return
      end
c-----------------------------------------------------------------------------
      SUBROUTINE XMAD (XDATA,N,XMED,XMADM)
c-----------------------------------------------------------------------------
c
c     The XMAD subroutine calculates the Median Absolute Deviation from
c     the sample median. The median, M , is subtracted from each
c     ORDERED statistic and then the absolute value is taken. This new
c     set of of statistics is then resorted so that they are ORDERED
c     statistics. The MAD is then defined to be the median of this
c     new set of statistics and is returned as XMADM. The MAD can
c     be defined:
c
c                   XMADM = median{ abs(x(i) - M) }
c
c     where the x(i) are the values passed in the array XDATA, and
c     the median, M, is passed in the array XLETTER. The set of stats
c     in the brackets is assumed to be resorted. For more information
c     see page 408 in UREDA.
c
c****************************************************************************
 
      parameter(iii=10000)
      implicit real*4 (a-h,o-z)
      dimension xdata2(iii),xdata(n)
      data dhalf,n1,n2/0.5,1,2/
      do i = 1,n
        xdata2(i) = abs(xdata(i) - xmed)
      enddo
      call qcksrt(n,xdata2)
ccc      call sort(n,xdata2)
      if (float(n)/float(n2) - int(n/n2) .eq. 0) then
        i1 = n/n2
        i2 = n/n2 + n1
        xmadm = dhalf*(xdata2(i1) + xdata2(i2))
      else
        i1 = int(n/n2) + n1
        xmadm = xdata2(i1)
      endif
      return
      end
c------------------------------------------------------------------------------
      SUBROUTINE MDIAN1 (X,N,XMED)
c------------------------------------------------------------------------------
c
c     Taken from Numerical Recipes, page 460.
c     Given an array X of N numbers, returns their median value XMED, and
c     the array is sorted low to high
c
c*******************************************************************************
      real x(n)
      call qcksrt(n,x)
ccc      call sort(n,x)
      n2=n/2
      if (2*n2.eq.n) then
        xmed=0.5*(x(n2)+x(n2+1))
      else
        xmed=x(n2+1)
      endif
      return
      end
*sbr mfs
*
*     computes the median, hinges, and fourth-spread of n observations
*
*     inputs:
*
*     x    =  input series, returned in ranked form
*     n    =  number of observations in series x
*
*     outputs:
*
*     xmed =  median of series x
*     hh   =  high hinge
*     hl   =  low hinge
*     fsd  =  fourth-spread
*     rmn  =  minimum value
*     rmx  =  maximum value
*
*     subroutines required:  sort
*
      subroutine mfs(x,n,xmed,hh,hl,fsd,rmn,rmx)
      real*4 x(n)
      integer temp1,temp2,j,k   !pjk
*
* rank the observations
      call qcksrt(n,x)
      rmn=x(1)
      rmx=x(n)
      k=n
      j=(k/2)+1
      temp1=n+1-j
      xmed=(x(j)+x(temp1))/2.0
*
* calculate the high (hh) and low (hl) hinges 
* for computing the fourth-spread (fsd)
      k=(k+1)/2
      j=(k/2)+1
      temp1=k+1-j
      hl=(x(j)+x(temp1))/2.0
      temp1=n-k+j
      temp2=n+1-j
      hh=(x(temp1)+x(temp2))/2.0
*
* the fourth-spread
      fsd=hh-hl
      return
      end
*sbr qcksrt
*
      subroutine qcksrt(n,arr)
      parameter (m=7,nstack=10000,fm=7875.,fa=211.,fc=1663.
     *    ,fmi=1.2698413e-4)
      dimension arr(n),istack(nstack)
      jstack=0
      l=1
      ir=n
      fx=0.
   10 if(ir-l.lt.m)then
        do 13 j=l+1,ir
          a=arr(j)
          do 11 i=j-1,1,-1
            if(arr(i).le.a)go to 12
            arr(i+1)=arr(i)
   11     continue
          i=0
   12     arr(i+1)=a
   13   continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        i=l
        j=ir
        fx=mod(fx*fa+fc,fm)
        iq=l+(ir-l+1)*(fx*fmi)
        a=arr(iq)
        arr(iq)=arr(l)
   20   continue
   21     if(j.gt.0)then
            if(a.lt.arr(j))then
              j=j-1
              go to 21
            endif
          endif
          if(j.le.i)then
            arr(i)=a
            go to 30
          endif
          arr(i)=arr(j)
          i=i+1
   22     if(i.le.n)then
            if(a.gt.arr(i))then
              i=i+1
              go to 22
            endif
          endif
          if(j.le.i)then
            arr(j)=a
            i=j
            go to 30
          endif
          arr(j)=arr(i)
          j=j-1
        go to 20
   30   jstack=jstack+2
        if(jstack.gt.nstack)pause 'nstack must be made larger.'
        if(ir-i.ge.i-l)then
          istack(jstack)=ir
          istack(jstack-1)=i+1
          ir=i-1
        else
          istack(jstack)=i-1
          istack(jstack-1)=l
          l=i+1
        endif
      endif
      go to 10
      end
*sbr cenbyseg - censor by minimum segment length - mnsl
*
      subroutine cenbyseg(mnsl,x,nc,idate,ident,mxy,mxs)
      real*4 x(mxy,mxs)
      integer idate(mxs,3)
      character ident(mxs)*8
      maxn=idate(nc+1,2)-idate(nc+1,1)+1
      nc1=0
      ifmny1=9999
      ilmxy1=-9999
c
c *** find those series .ge. minimum segment legnth and save
c
      do j=1,nc
        nyrs=idate(j,3)
        if(nyrs.ge.mnsl)then
          nc1=nc1+1
          idate(nc1,1)=idate(j,1)
          idate(nc1,2)=idate(j,2)
          idate(nc1,3)=idate(j,3)
          ident(nc1)=ident(j)
          if(idate(nc1,1).lt.ifmny1)ifmny1=idate(nc1,1)
          if(idate(nc1,2).gt.ilmxy1)ilmxy1=idate(nc1,2)
          do i=1,maxn
            x(i,nc1)=x(i,j)
          enddo
        endif
      enddo
      maxn1=ilmxy1-ifmny1+1
      idate(nc1+1,1)=ifmny1
      idate(nc1+1,2)=ilmxy1
      idate(nc1+1,3)=maxn1
      if(maxn1.lt.maxn)then
        ndif=maxn-maxn1
        write(*,'(/1x,''Total chronology length shortened: '',i6,''  years'')')ndif
c        do i=1,maxn1
c          do j=1,nc
c            x(i,j)=x(i+ndif,j)
c          enddo
c        enddo
      endif
c
c *** dummy out the data matrix with missing values to account for removed series
c
      do j=nc1+1,nc
        do i=1,maxn
          x(i,j)=-9.99
        enddo
        call blank(ident(j))
      enddo
c
c *** final results
c
      write(*,'(/t30,''    NS   IFY   ILY   NYR'')')
      write(*,'(1x,''Min segment length series:'',t30,4i6,''  Years total time span'')')
     *  nc1,ifmny1,ilmxy1,maxn1
      nc=nc1
      return
      end
*sbr cenbyify - censor by maximum starting year - mxsy
*
      subroutine cenbyify(mxsy,x,nc,idate,ident,mxy,mxs)
      real*4 x(mxy,mxs)
      integer idate(mxs,3)
      character ident(mxs)*8
      maxn=idate(nc+1,2)-idate(nc+1,1)+1
      nc1=0
      ifmny1=9999
      ilmxy1=-9999
c
c *** find those series .le. minimum starting year and save
c
      do j=1,nc
        ify=idate(j,1)
        if(ify.le.mxsy)then
          nc1=nc1+1
          idate(nc1,1)=idate(j,1)
          idate(nc1,2)=idate(j,2)
          idate(nc1,3)=idate(j,3)
          ident(nc1)=ident(j)
          if(idate(nc1,1).lt.ifmny1)ifmny1=idate(nc1,1)
          if(idate(nc1,2).gt.ilmxy1)ilmxy1=idate(nc1,2)
          do i=1,maxn
            x(i,nc1)=x(i,j)
          enddo
        endif
      enddo
      maxn1=ilmxy1-ifmny1+1
      idate(nc1+1,1)=ifmny1
      idate(nc1+1,2)=ilmxy1
      idate(nc1+1,3)=maxn1
      if(maxn1.lt.maxn)then
        ndif=maxn-maxn1
        write(*,'(/1x,''Total chronology length shortened: '',i6,''  years'')')ndif
c        do i=1,maxn1
c          do j=1,nc
c            x(i,j)=x(i+ndif,j)
c          enddo
c        enddo
      endif
c
c *** dummy out the data matrix with missing values to account for removed series
c
      do j=nc1+1,nc
        do i=1,maxn
          x(i,j)=-9.99
        enddo
        call blank(ident(j))
      enddo
c
c *** final results
c
      write(*,'(/t30,''    NS   IFY   ILY   NYR'')')
      write(*,'(1x,''Max starting year series:'',t30,4i6,''  Years total time span'')')
     *  nc1,ifmny1,ilmxy1,maxn1
      nc=nc1
      return
      end
*sbr cenbycom - censor by common period - ifc, ilc
*
      subroutine cenbycom(ifc,ilc,x,nc,idate,ident,mxy,mxs)
      real*4 x(mxy,mxs)
      integer idate(mxs,3),ifc,ilc,nc
      character ident(mxs)*8
      maxn=idate(nc+1,2)-idate(nc+1,1)+1
      nc1=0
      nyc=ilc-ifc+1
      ifmny1=ifc
      ilmxy1=ilc
      maxn1=ilmxy1-ifmny1+1
c
c *** find those series with data in the common period - ifc, ilc
c
      do j=1,nc
        ify=idate(j,1)
        ily=idate(j,2)
        if(ify.le.ifc.and.ily.ge.ilc)then
          nc1=nc1+1
          idate(nc1,1)=idate(j,1)
          idate(nc1,2)=idate(j,2)
          idate(nc1,3)=idate(j,2)-idate(j,1)+1
          ident(nc1)=ident(j)
          do i=1,maxn
            x(i,nc1)=x(i,j)
          enddo
        endif
      enddo
      idate(nc1+1,1)=ifc
      idate(nc1+1,2)=ilc
      idate(nc1+1,3)=nyc
c
c *** shorten the length of retained series to the common interval only
c
      write(*,'(/1x,''Total chronology length shortened: '',i6,''  years'')')maxn-nyc
      do j=1,nc1
        ify=idate(j,1)
        ily=idate(j,2)
        ndif=ifc-ify
        do i=1,nyc
          x(i,j)=x(i+ndif,j)
        enddo
        do i=1,ndif
          x(i+nyc,j)=-9.99
        enddo
        idate(j,1)=ifc
        idate(j,2)=ilc
        idate(j,3)=nyc
      enddo
c
c *** dummy out the data matrix with missing values to account for removed series
c
      do j=nc1+1,nc
        do i=1,maxn
          x(i,j)=-9.99
        enddo
        call blank(ident(j))
      enddo
c
c *** final results
c
      write(*,'(/t30,''    NS   IFY   ILY   NYR'')')
      write(*,'(1x,''Max starting year series:'',t30,4i6,''  Years total time span'')')
     *  nc1,ifc,ilc,nyc
      nc=nc1
      return
      end
*sbr savcen - save the censored data in Tucson raw ringwidth format
*
      subroutine savcen(cproc,icen1,icen2,x,y,nc,idate,ident,itrdb,filnam,iu11,mxy,mxs)
      real*4 x(mxy,mxs),y(mxy)
      integer idate(mxs,3),icen1,icen2,nc
      character cproc,filnam*40,fil1*50,ident(mxs)*8,id*8,itrdb(3)*85,tag*3
c
c *** parse the input file name for creating the output file name
c
      ll=len(filnam)
      jj=0
      do j=1,ll
        jj=ll-j+1
        if(filnam(jj:jj).ne.' ')goto 1
      enddo
    1 ll=jj
      do j=1,ll
        jj=j
        if(filnam(j:j).eq.'.')then
          jj=j-1
          goto 2
        endif
      enddo
    2 ll=jj
c
c *** create the unique output file name for the censored data type
c
      call blank(fil1)
      fil1(1:40)=filnam(1:40)
      if(cproc.eq.'2')then
        if(icen1.lt.100)then
          fil1(ll+1:ll+12)='_mnsl=  .txt'
          write(fil1(ll+7:ll+8),'(i2)')icen1
        elseif(icen1.ge.100.and.icen1.lt.1000)then
          fil1(ll+1:ll+13)='_mnsl=   .txt'
          write(fil1(ll+7:ll+9),'(i3)')icen1
        elseif(icen1.ge.1000)then
          fil1(ll+1:ll+14)='_mnsl=    .txt'
          write(fil1(ll+7:ll+10),'(i4)')icen1
        endif
      elseif(cproc.eq.'3')then
        fil1(ll+1:ll+14)='_mxsy=    .txt'
        write(fil1(ll+7:ll+10),'(i4)')icen1
      elseif(cproc.eq.'4')then
        fil1(ll+1:ll+19)='_comm=    -    .txt'
        write(fil1(ll+7:ll+10),'(i4)')icen1
        write(fil1(ll+12:ll+15),'(i4)')icen2
      endif
      call blank(tag)
      if(cproc.eq.'2')tag='seg'
      if(cproc.eq.'3')tag='syr'
      if(cproc.eq.'4')tag='com'
      open(iu11,file=fil1,status='unknown')
      write(*,'(/1x,''File for censored data created: '',a)')fil1
      do j=1,3
        if(cproc.eq.'2')itrdb(j)(81:85)=' MNSL'
        if(cproc.eq.'3')itrdb(j)(81:85)=' MXSY'
        if(cproc.eq.'4')itrdb(j)(81:85)=' COMM'
        write(iu11,'(a)')itrdb(j)
      enddo
c
c *** write the censored data to the file in Tucson ring width format
c
      do j=1,nc
        call blank(id)
        id=ident(j)
        ify=idate(j,1)
        nyr=idate(j,3)
        do i=1,nyr
          y(i)=x(i,j)
        enddo
        call trrw(id,ify,nyr,y,nss,ifl,iu11,iuo,'mw',tit,tag,mxy,itrdb)
      enddo
      close(iu11)
      return
      end