*
      Program RCSsigFree
c
c *** does the minimum and implements Tom Melvin's signal free method
c     and even more minimal RCS detrending
c
c Subroutine list
c      subroutine sigfree(x,y,f,ls,alpha,sfind,nc,crn1,crn2,crn3,nss,idate,ident,iopt,irob,
c     *                     switch,reset,hpfilt,tol,mxy,mxs,itrdb,ybar,rxm,runr,nsg,filnam,fil1,rcs,pith)
c      subroutine curvefits(y,f,iyr,n,n1,id,iopt,alpha,mxy,mxs)
c      subroutine curve(n,n1,z,f,eq)		       !hakan
c      subroutine trend(n,n1,y,f,sl,yi)			!hakan
c      subroutine xmean(x,f,n,n1)
c      subroutine supsmu (n,x,y,w,iper,span,alpha,smo,sc)
c      subroutine smooth (n,x,y,w,span,iper,vsmlsq,smo,acvr)
c      subroutine spline(n,y,f,pp,ls)
c      subroutine spline1(n,y,f,ls)
c      subroutine spline2(n,y,f,ls)
c      subroutine robsum(x,n,nc,crnmn,crnsd,nss,iopt,mxy,mxs)
c      subroutine mnsd(x,n,xm,sd)
c      subroutine xbiwt (xdata,n,xlbiwt,xsbiwt,xlbiwt1,xsbiwt1)
c      subroutine xmad (xdata,n,xmed,xmadm)
c      subroutine mdian1 (x,n,xmed)
c      subroutine mfs(x,n,xmed,hh,hl,fsd,rmn,rmx)
c      subroutine qcksrt(n,arr)
c      subroutine cenbyseg(mnsl,x,nc,idate,ident,mxy,mxs)
c      subroutine cenbyify(mxsy,x,nc,idate,ident,mxy,mxs)
c      subroutine cenbycom(ifc,ilc,x,nc,idate,ident,mxy,mxs)
c      subroutine savcen(cproc,icen1,icen2,x,y,nc,idate,ident,itrdb,filnam,iu11,mxy,mxs)
c      subroutine arpad(x,n,y,np,ls,mxy)
c      subroutine memcof(data,n,m,pm,cof,wk1,wk2,wkm)
c      subroutine extend(y,n,ncast,arps,ip,iopt)
c      subroutine plot1(nmax,wrk,sp,ny,iy,xlab,ylab,tit,rline,annot1,annot2)
c      subroutine spaghetti(xy,ny,ns,ify,xlab,ylab,tit,annot1,annot2,mxy,mxs)
c      subroutine runbar(x,ybar,idate,nc,rxm,runr,ncnt,nsg,nlp,mxy,mxs,iu33)
c      subroutine rrbar(x,wrk1,idate,nc,ifskp,ilskp,nr,ave,sdev,see,nsg,mxy,mxs,maxlen)
c      subroutine pearsn(x,y,n,r,a,b)
c      subroutine intfrm(datum,frm)
c      subroutine stringclean(string,size)
c      subroutine trrw(id,iyr,n,y,nc,ifl,iu,iuo,q,tit,tag,mxy,itrdb)
c      subroutine datar(id,n,x,iyr,iu,mxy)
c      subroutine blank(l)
c      subroutine trread(x,y,nc,idate,ident,itrdb,iu10,mxy,mxs)
c RCS subroutines and calls
c      subroutine plt1(mxsl,rmax,mey,ycal,fycaly,lycaly,flag)
c	subroutine plt2(mxy,mxs,pltdata,nc,idate,flag)
c	subroutine plt3(mxy,mxs,indx,xrc,f,idate,flag,nmax,sernum,n,ident)
c	subroutine plt4(mxlen,maxit,rcscrv,nmax,rcnt)
c      subroutine rcurve(mxy,mxs,n,x,xrc,f,idate,pith,flag,mxage,eq,ndt)
c      subroutine pithy(mxy,iu,pith,poff)
c      subroutine detrend2(x,f,n1,n,eq,jdt,mxy,ndt)			!Hakan
c      subroutine gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n1,ne,x,y,yfit,sumsq)	!Hakan
c      subroutine amoeba(p,y,mp,np,ndim,iter,ftol,iopt,
c      subroutine hughdi(x,y,n1,ntot,xy,coef,flag)		!Hakan
c      subroutine gexp(y,nt,n,f,wk1,wk2,coef,flag)		!Hakan
c      subroutine moment(data,n,ave,adev,sdev,var,skew,curt,sens,ac)
c      subroutine xmean2(x,n,xm)
c      subroutine lowess(mxd,y,ys,n,f)   ! added mxd to argument list
c      subroutine lowest(x,y,n,xs,ys,nleft,nright,w,userw,rw,ok)
c      subroutine msmoor(x,y,z,n,m,iter,ier,mxy)
c      subroutine msmoo(x,y,z,n,m,iter,ier,mxy)
c      subroutine sort(n,ra)
c      function gamoff(v,ndmax,nemax,ne,x,y,iopt)
c      function ametry(p,y,psum,mp,np,ndim,ihi,fac,iopt,nemax,ne,xe,ye)
c      function minimum(inarray,arraynum)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      parameter(mxy=5000,mxs=3000)
      parameter(iu10=10,iu11=11,iu12=12,iu69=69)
      real*4 x(mxy,mxs),y(mxy),f(mxy),sfind(mxy,mxs)
      real*4 crn1(mxy,2),crn2(mxy,2),crn3(mxy,2)
      real*4 ybar(mxy,3),rxm(mxy,7),runr(mxy,5)
      real*8 tol
      integer iopt,irob
      integer idate(mxs,3),nss(mxy),ls(mxy)
      character filnam*40,fil1*60,itrdb(3)*88,ident(mxs)*8,id*8
ccc      character tit*85,xlab*20,ylab*20,annot1*35,annot2*35
      character itrun,melvin,switch,reset,hpfilt,itol,cproc,cdat*10,frm*20
c new globasl decs pjk
	integer ndt
	integer pith(mxy)
	logical poff,rcs,tryagain,fok
c		
	rcs=.false.;poff=.false.;tryagain=.false.;fok=.false. 
c
c *** open data file to process
c
      write(*,'(/t24,''**********************************''/
     *           t24,''**********************************''/
     *           t24,''***                            ***''/
     *           t24,''***     Program RCSsigFree     ***''/
     *           t24,''***                            ***''/
     *           t24,''***       Creation Date:       ***''/
     *           t24,''***         05/11/2011         ***''/
     *           t24,''***                            ***''/
     *           t24,''***       Programmed by:       ***''/
     *           t24,''***  Edward R. Cook: TRL-LDEO  ***''/
     *           t24,''***  Paul J. Krusic: BBCC-SU   ***''/
     *           t24,''***  Tom Melvin: CRU-UEA       ***''/
     *           t24,''***                            ***''/
     *           t24,''********** Version 2.2a **********''/
     *           t24,''**********************************'',/)')
c
112   call blank(filnam)
      write(*,'(/1x,''Tree-ring data file name (/ to exit)'',t40,''==> '',$)')
      read(*,'(a)')filnam
      if(filnam(1:1).eq.'/')stop
      inquire(file=filnam,exist=fok)
      if(.not.fok)then
	  write(*,'(//10x,''File > '',a25,'' < not found'')')filnam
	  pause'         Hit return to enter a new file name or quit program.'
	  goto 112
      endif
      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)=-999.0
          sfind(i,j)=-999.0
        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
   10 call blank(cdat)
      write(*,'(//1x,''General Curve fit options:''//
     *  1x,''Spline pct N filter (variable)  .le. -10''/
     *  1x,''Friedman variable span smoother'',t40,''-5''/
     *  1x,''Regional Curve Standardization (RCS)'',t40,''-2''/
     *  1x,''Neg Expon/linear any slope'',t40,''<1>''/
     *  1x,''Neg Expon/linear neg slope'',t41,''2''/
     *  1x,''Linear regression any slope'',t41,''3''/
     *  1x,''Linear regression neg slope'',t41,''4''/
     *  1x,''Horizontal line through mean'',t41,''5''/
c     *  1x,''Age-dependent spline smoothing'',t41,''8''/
     *  1x,''Median segment length spline'',t41,''9''/
     *  1x,''Spline filter in years (fixed)  .ge.  10''/
     *  1x,''Opt 1 is default, hit <RET> if you want it'',t50,''==> '',$)')
      read(*,'(a)')cdat
      if(cdat(1:1).eq.'/')stop
      if(cdat(1:1).eq.' ')cdat(1:1)='1'
      call intfrm(cdat,frm)
      read(cdat,frm,err=10)iopt
      if(iopt.eq.9)iopt=medseg
      if(iopt.eq.-5)then
        write(*,'(/1x,''Enter Friedman tweeter (alpha) 0-9'',t50,''==> '',$)')
        read(*,*,err=10)alpha
        if(alpha.lt.0.or.alpha.gt.9)then
          write(*,'(/1x,''Friedman tweeter not in admissible range'')')
          goto 10
        endif
      endif
c -----------------------------------------------------------------------------new pjk
      if(iopt.eq.-2)then
          rcs=.true.;tryagain=.false.
111       write(*,'(//1x,''Select RCS curve fit options:'',//)')
          write(*,'(t3,''  -5:  friedman super smoother'',
     *      ''        lower alpha (0-9) gives more local fit'')')
          write(*,'(t3,''  -4:  robust lowess smoother'',
     *      ''         integer percent of series length <Not Yet>'')')
          write(*,'(t3,''  -3:  median smoothing'',
     *      ''               must have odd window length <full length only>'')')
          write(*,'(t3,''  -2:  Age-dependent spline smoothing'')')
          write(*,'(t3,''   1:  neg expon curve, '',
     *      ''no = opt 4    f(i) = a*exp(-b*t(i)) + d    (d>0)'')')
          write(*,'(t3,''   2:  neg expon curve, '',
     *      ''no = opt 5    f(i) = a*exp(-b*t(i)) + d    (d>0)'')')
          write(*,'(t3,''   3:  neg expon curve  '',
     *      ''(general)     f(i) = a*exp(-b*t(i+p)) + d (any d)'')')
          write(*,'(t3,''   4:  linear regression '',
     *      ''(any slope)  f(i) = +/-c*t(i) + d'')')
          write(*,'(t3,''   5:  linear regression '',
     *      ''(neg slope)  f(i) = -c*t(i) + d   '')')
          write(*,'(t3,''   6:  horizontal line through '',
     *      ''mean   f(i) = mean(y(i)) = d'')')
          write(*,'(t3,''   7:  hugershoff growth '',
     *      ''function     f(i) = a*t(i+p)**b * exp(-c*t(i+p)) + d'')')
          write(*,'(t3,''   8:  general exponential (b = 1)    '',
     *      ''f(i) = a*t(i+p)**b * exp(-c*t(i+p)) + d'')')
          write(*,'(t3,''  >9: cubic smoothing '',
     *      ''spline         fixed 50 pct variance cutoff'')')
          write(*,'(t3,'' <-9: cubic smoothing '',
     *      ''spline         pct n 50 pct variance cutoff'')')
          write(*,'(/1x,''enter  option  (x to start over)  ==> '',$)')
          read(*,*,err= 10)iopt
	   if(iopt.eq.-5)then
	        ndt=0
	        write(*,'(/1x,''Enter tweeter alpha 0-9'',t50,''==> '',$)')
	        read(*,*)ndt
	        if(ndt.lt.0.or.ndt.gt.9)then
			write(*,'(/1x,''Friedman alpha not in admissible 0-9 range'')')
			goto 111
		endif
	   elseif(iopt.eq.-4)then
              ndt=0
              write(*,'(/1x,''enter the lowess window width''/
     *          1x,''fixed for all series: a positive integer >2''/
     *          1x,''and less then the shortest series length, or''/
     *          1x,''percent of series length: a negative integer''/
     *          1x,''between -1 and -100, inclusive'', t45,''==> '',$)')
              read(*,*)ndt
              if(ndt.ge.0)then
                  if(ndt.eq.0.or.ndt.eq.1.or.ndt.eq.2)goto 111
              else
                  if(ndt.gt.-1)goto 111
                  if(ndt.lt.-100)goto 111
              endif
	   elseif(iopt.eq.-3)then
              ndt=0
              write(*,'(/1x,''enter the median window width''/
     *          1x,''fixed for all series: a positive odd integer >2''/
     *          1x,''and less then the shortest series length, or''/
     *          1x,''percent of series length: a negative integer''/
     *          1x,''between -1 and -100, inclusive'', t45,''==> '',$)')
              read(*,*)ndt
              if(ndt.ge.0)then
          	    if(ndt.eq.0.or.ndt.eq.1.or.ndt.eq.2)goto 111
          	    if(mod(ndt,2).eq.0)goto 111
        	else
          	    if(ndt.gt.-1)goto 111
          	    if(ndt.lt.-100)goto 111
        	endif
	   endif
	   alpha=real(ndt)
	   if(.not.tryagain)call pithy(mxs,iu69,pith,poff)
      endif  !iopt =-2
c ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^new pjk
c
c *** calculate mean chronology robustly or stupidly not
c
   11 call blank(cdat)
      write(*,'(/1x,''How to calculate the mean chronology (/ to stop):''/
     *           1x,''Non-robust arithmetic mean'',t39,'' 1 ''/
     *           1x,''Robust Tukey biweight mean'',t39,''<2>''/
     *           1x,''Robust is default, hit <RET> if you want it'',t50,''==> '',$)')
      read(*,'(a)')cdat
      if(cdat(1:1).eq.'/')stop
      if(cdat(1:1).eq.' ')cdat(1:1)='2'
      if(cdat(1:1).eq.'1')then
        irob=1
      elseif(cdat(1:1).eq.'2')then
        irob=2
      else
        write(*,'(/1x,''Mean index choice not admissible'')')
        goto 11
      endif
c
c *** get running RBAR window width for performing running RBAR variance stabilization
c
   12 call blank(cdat)
      write(*,'(/1x,''Enter running RBAR window width for variance''/
     *           1x,''stabilization. Must be .ge. 25 years and odd.''/
     *           1x,''Default window width is 51 years, 0 to skip'',t50,''==> '',$)')
      read(*,'(a)')cdat
      if(cdat(1:1).eq.'/')stop
      if(cdat(1:1).eq.' ')cdat(1:2)='51'
      call intfrm(cdat,frm)
      read(cdat,frm,err=12)nsg
      if(nsg.eq.0)then
        write(*,'(/1x,''Running RBAR variance stabilization skipped'')')
      elseif(nsg.gt.0.and.nsg.lt.25)then
        write(*,'(/1x,''Running RBAR window width must be .ge. 25'')')
        goto 12
      elseif(mod(nsg,2).eq.0)then
        nsg=nsg+1
        write(*,'(/1x,''Running RBAR window width increased to make it odd: '',i5)')nsg
      endif
c
      itrun='n'
c
   13 melvin=' '
      switch='n'
      reset='y'
      hpfilt='y'
      tol=0.0001d0
      write(*,'(/1x,''Play with the signal free convergence settings?''/
     *           1x,''WARNING: This is a very dangerous thing to do!''/
     *           1x,''Enter y/<n> (default is n for NO, / to exit)'',t50,''==> '',$)')
      read(*,'(a)')melvin
      if(melvin.eq.'/')stop
      if(melvin.eq.' ')melvin='n'
      if(melvin.eq.'N')melvin='n'
      if(melvin.eq.'Y')melvin='y'
      if(melvin.ne.'y'.and.melvin.ne.'n')then
        write(*,'(/1x,''Melvin not entered as y or n'')')
        goto 13
      endif
      
      if(melvin.eq.'y')then
        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
        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
        hpfilt=' '
        write(*,'(/1x,''H-P filter chronology differences? <y>/n'',t50,''==> '',$)')
        read(*,'(a)')hpfilt
        if(hpfilt.eq.'/')stop
        if(hpfilt.eq.' ')hpfilt='y'
        if(hpfilt.eq.'Y')hpfilt='y'
        if(hpfilt.eq.'N')hpfilt='n'
   20   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:   ''/
     *             1x,''  Original Melvin setting   1  0.0005''/
     *             1x,''  Current default setting  <2> 0.0001''/
     *             1x,''  **** Roll Your Own ****   3  X.XXXX'',t50,''==> '',$)')
        read(*,'(a)')itol
        if(itol.eq.'/')stop
        if(itol.eq.'1')then
          tol=0.0005d0
        elseif(itol.eq.' '.or.itol.eq.'2')then
          tol=0.0001d0
        elseif(itol.eq.'3')then
          write(*,'(/1x,''Okay, enter YOUR convergence test'',t50,''==> '',$)')
          read(*,*,err=20)tol
        else
          goto 20
        endif
      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,crn3,nss,idate,ident,iopt,irob,
     *  switch,reset,hpfilt,tol,mxy,mxs,itrdb,ybar,rxm,runr,nsg,filnam,fil1,rcs,pith,ndt,tryagain)
     
      if(tryagain)goto 111
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)(83:88)='stdcrn'
        write(iu11,'(a)')itrdb(j)
      enddo
      call blank(id)
      id=filnam(1:6)
      call trrw(id,ifmny,maxn,crn1(1,1),nss,ifl,iu11,iuo,'iw',tit,'stdcrn',mxy,itrdb)
      if(nsg.gt.0)then
        do j=1,3
          itrdb(j)(83:88)='stdstb'
          write(iu11,'(a)')itrdb(j)
        enddo
        call blank(id)
        id=filnam(1:6)
        call trrw(id,ifmny,maxn,crn3(1,1),nss,ifl,iu11,iuo,'iw',tit,'stdstb',mxy,itrdb)
      endif
      do j=1,3
        itrdb(j)(83:88)='ssfcrn'
        write(iu11,'(a)')itrdb(j)
      enddo
      call trrw(id,ifmny,maxn,crn2(1,1),nss,ifl,iu11,iuo,'iw',tit,'ssfcrn',mxy,itrdb)
      if(nsg.gt.0)then
        do j=1,3
          itrdb(j)(83:88)='ssfstb'
          write(iu11,'(a)')itrdb(j)
        enddo
        call blank(id)
        id=filnam(1:6)
        call trrw(id,ifmny,maxn,crn3(1,2),nss,ifl,iu11,iuo,'iw',tit,'ssfstb',mxy,itrdb)
      endif
      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')
      if(nsg.gt.0)then
        write(iu11,'(''  Year  SSize  StdCrn  SsfCrn  StdSig  SsfSig  StdSee  SsfSee  StdStb  SsfStb'')')
      else
        write(iu11,'(''  Year  SSize  StdCrn  SsfCrn  StdSig  SsfSig  StdSee  SsfSee'')')
      endif
      do i=1,maxn
        rns=float(nss(i))
        if(nsg.gt.0)then
          write(iu11,'(2i6,1x,12f8.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),crn3(i,1),crn3(i,2)
        else
          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)
        endif
      enddo
      close(iu11)
  900 write(*,'(/1x,''Normal Termination of program'',/1x,''Hit RETURN to Quit'',$)')
  	pause
      stop
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr sigfree
*
      subroutine sigfree(x,y,f,ls,alpha,sfind,nc,crn1,crn2,crn3,nss,idate,ident,iopt,irob,
     *  switch,reset,hpfilt,tol,mxy,mxs,itrdb,ybar,rxm,runr,nsg,filnam,fil1,rcs,pith,ndt,tryagain)
     
c
      parameter(mxlen=5000,mxser=3000,maxit=100,iu29=29,iu30=30,iu31=31,iu32=32,iu33=33)
      real*4 x(mxy,mxs),y(mxy),f(mxy),sfind(mxy,mxs),crn1(mxy,2),crn2(mxy,2),crn3(mxy,2)
      real*4 crnsav(mxlen),xbar(mxser),segmn(mxlen),segsd(mxlen),agemn(mxlen),agesd(mxlen)
      real*4 test1(maxit),test2(maxit),test3(maxit),test4(maxit),wrk1(mxlen),wrk2(mxlen)
      real*4 ybar(mxy,3),rxm(mxy,7),runr(mxy,5),savit(mxlen,maxit),ones(mxlen)
      real*8 test,tol,testd,testdd      
      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,hpfilt,line*410,reset,switch
      character star1,star2,star3,star4,itrdb(3)*88,id*8,istab,filnam*40,fil1*60
c new RCS declarations
      integer p,j,k,jj
      real*4 xrc(mxlen),xtmp(mxlen,mxser),eq(7),yrcs(mxlen,mxser),rcscrv(mxlen,2,maxit),flag,miny    !sample regional curve
      parameter (flag=-999.0)
      integer pith(mxy),mxage,ndt,rcmax,rcnt,age(mxy),fyr(mxser)
      logical rcs,tryagain,rcscut,itype,plotall,pt
      character ans*1, string*35
c
      if(nsg.eq.0)then
        istab='n'
      else
        istab='y'
      endif
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
        age(j)=n              !pjk
        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:37)='Mean Segment Length And Mean Tree Age'
      call blank(xlab)
      xlab(1:5)='Years'
      call blank(ylab)
      ylab(1:5)='Years'
      call blank(annot1)
      annot1(1:20)=' Mean Segment Length'
      call blank(annot2)
      annot2(1:20)='Mean Approx Tree 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
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
      call blank(fil1)
      fil1(1:40)=filnam(1:40)
      fil1(ll+1:ll+15)='_OrigCurves.txt'
      open(iu31,file=fil1,status='unknown')
      close(iu31,status='delete')
      open(iu31,file=fil1,status='unknown')
      fil1(ll+1:ll+16)='_OrigIndices.txt'
      open(iu32,file=fil1,status='unknown')
      close(iu32,status='delete')
      open(iu32,file=fil1,status='unknown') 
c___________________________________________ RCS       
      itype=.false.; pt=.false.
111   write(*,'(//1x,''Use [r]atios or residual[s]? : '',$)')
      read(*,'(a1)',err=111)ans
      if((ans.eq.'S').or.(ans.eq.'s'))then
         itype=.true.
1111     write(*,'(//1x,''Power Transform data? [t]rue or [f]alse : '',$)')
         read(*,*,err=1111)pt
         if(pt)then  
           do k=1,nc
             p=1
             do j=idate(k,1),idate(k,2)
               xtmp(p,k)=x(p,k) !copy the original input data matrix
               p=p+1
             enddo
           enddo
           call slpt(mxy,mxs,x,idate,nc) !pt
         endif
      endif
c
      if(rcs)then  ! new  rcs
         rcmax=0
112      call rcurve(5000,3000,nc,x,xrc,f,idate,pith,flag,mxage,eq,iopt,alpha,ndt,poff,rcmax) !(xrc is the regional mean (f is the regional curve)
	  call plt1(mxy,mxs,xrc,f,1,mxage,flag,'Initial RC and MVF fit             ')
c        save the original Regional Curve and it's solution
	  rcnt=1
	  do jj=1,mxage
	     rcscrv(jj,1,rcnt)=xrc(jj);rcscrv(jj,2,rcnt)=f(jj)
	  enddo
	  tryagain=.false.
	  do j=1,mxage
	     if(f(j).le.0.and.f(j).ne.flag)then
	        tryagain=.true.
	        goto 1113
	    endif
	  enddo
1113	  if(tryagain)then
	        write(*,'(/3x,''**************** WARNING ****************'')')
	        write(*,'(3x,''    Regional Curve goes south (negative) '')')
	        write(*,'(3x,'' Suggest choosing another RCS configuration.'')')
	        write(*,'(/3x,''**************** WARNING ****************'')')
	        pause '                   Hit RETURN to Continue'
	  endif
	  write(*,'(//1x,''Try a different RCS detrending method?'')')
	  write(*,'(1x,'' [t]rue or [f]alse : '',$)')
	  read(*,*)tryagain
	  if(tryagain.and.pt)then
             do k=1,nc
                 p=1
                 do j=idate(k,1),idate(k,2)
                   x(p,k)=xtmp(p,k)
                   p=p+1
                 enddo
             enddo
	      return
	  elseif(tryagain.and.(.not.pt))then
	      return
	  endif
	  write(*,'(//1x,''Truncate Curve?'')')
	  write(*,'(1x,'' [t]rue or [f]alse : '',$)')
	  read(*,*)rcscut
	  if(rcscut)then
	     do j=1,nc
	        fyr(j)=idate(j,1)
	     enddo
	     call comint(mxser,fyr,ident,age,nc)
	     write(*,'(//3x,''Enter Biological Year beyond which not to use'')')
	     write(*,'(3x,''in the estimation of the Regional model : '',$)')
	     read(*,*)rcmax
	     goto 112
	  endif
	  write(*,'(/1x,''Make lots of plots [t]rue or [f]alse ? : '',$)')
	  read(*,*)plotall
      endif        ! new  rcs
      miny=9999.
c___________________________________________  end RCS      
      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
c
          if(.not.rcs)then                    
             call curvefits(y,f,ify,n,n1,tit,iopt1(j,2),alpha,mxy,mxs)
          endif             
          do i=1,n   !compute indices or residuals from series y
             if(.not.itype)then
                wrk2(i)=y(i)/f(i)
                sfind(i,j)=y(i)/f(i)   !indices
             else
                wrk2(i)=y(i)-f(i)
                sfind(i,j)=y(i)-f(i)   !residuals
                if(wrk2(i).lt.miny)miny=wrk2(i)          !min resid
             endif
          enddo
          if(.not.itype)then
             call trrw(ident(j),ify,n,f,nss,ifl,iu31,iuo,'mw',tit,'stdcrv',mxy,itrdb)
             call trrw(ident(j),ify,n,wrk2,nss,ifl,iu32,iuo,'mw',tit,'stdind',mxy,itrdb)
          endif
      enddo
c
      if(itype)then  !if residuals convert all to positive values by adding (abs(miny) + 0.1) to all residuals
          if(miny.ge.0.)then
             miny=0.000
          elseif(miny.lt.0)then
             miny=abs(miny)+0.00100    !scaling factor for negative residuals
          endif
          do j=1,nc
             n=idate(j,2)-idate(j,1)+1
             do i=1,n
                sfind(i,j)=sfind(i,j)+miny
                wrk2(i)=sfind(i,j)
             enddo
             call trrw(ident(j),ify,n,f,nss,ifl,iu31,iuo,'mw',tit,'stdcrv',mxy,itrdb)
             call trrw(ident(j),ify,n,wrk2,nss,ifl,iu32,iuo,'mw',tit,'stdind',mxy,itrdb)
          enddo
      endif
      if(rcs)then
         if(itype)then
            string='Initial RC Residual Indices        '
         else
            string='Initial RC Ratio Indices           '
         endif
         call plt2(mxy,mxs,sfind,nc,idate,flag,'Initial RC and MVF fit             ')
      endif
c
      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)=-999.0
            sfind(i,j)=-999.0
          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)
      call mnsd(crn1(1,1),maxn,xmstd1,sdstd1)
      do i=1,maxn
        ybar(i,1)=crn1(i,1)				!running correlation bit added 01/03/11
        runr(i,1)=float(nss(i))
      enddo
      nss1=0
      do i=1,maxn
        if(nss(i).eq.1)then
          nss1=nss1+1
        else
          goto 5
        endif
      enddo
    5 write(*,'(/1x,''nss1 skip first number of years: '',i3)')nss1
      if(nss1.gt.0)then
        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
      else
        nreset=maxn
      endif

      call mnsd(crn1(nss1+1,1),nreset,xm,sd)
      do i=1,maxn
ccc        crn1(i,1)=crn1(i,1)-xm+1.0
        crn2(i,1)=crn1(i,1)
        crnsav(i)=crn1(i,1)
      enddo
c
c *** perform running rbar variance stabilization on original standard indices
c
      if(istab.eq.'y')then
        write(*,'(/1x,''Perform standard indices running rbar calculations'')')
        nlp=nsg-1
        call blank(fil1)
        fil1(1:40)=filnam(1:40)
        fil1(ll+1:ll+13)='_Std_RBAR.txt'
        open(iu33,file=fil1,status='unknown')
        call runbar(sfind,ybar,idate,nc,rxm,runr,ncnt,nsg,nlp,mxy,mxs,iu33)
        close(iu33)
        call blank(tit)
        tit(1:45)='Standard Indices Running Average Correlations'
        call blank(xlab)
        xlab='Years'
        call blank(ylab)
        ylab='Correlations'
        call blank(annot1)
        annot1(1:14)='Running RBAR ('
        if(nsg.lt.100)then
          write(annot1(15:16),'(i2)')nsg
          annot1(17:21)=' Yrs)'
        else
          write(annot1(15:17),'(i3)')nsg
          annot1(18:22)=' Yrs)'
        endif
        call blank(annot2)
        annot2(1:20)='Constant Sample RBAR'
        ones(1)=-999.9
        call plot1(mxy,runr(1,2),runr(1,4),maxn,ifmny,xlab,ylab,tit,ones,annot1,annot2)
        do i=1,maxn
          crn3(i,1)=ybar(i,3)
        enddo
        call mnsd(crn3(1,1),maxn,xmstd2,sdstd2)
        do i=1,maxn
          datum=crn3(i,1)-xmstd2
          crn3(i,1)=datum+xmstd1
          ones=1.0
        enddo
        call blank(tit)
        tit(1:53)='Original vs. Running RBAR Stabilized Standard Indices'
        call blank(annot1)
        annot1(1:23)='Original STD Chronology'
        call blank(annot2)
        annot2(1:21)='Stable STD Chronology'
        call plot1(mxy,crn1(1,1),crn3(1,1),maxn,ifmny,xlab,ylab,tit,ones,annot1,annot2)
      endif
c
c *** do the signal free iterations
c
      iter=0
      small=9999.0
      smooth1=9999.0
      smooth2=9999.0
      big=9999.0
c      write(*,'(/38x,''Mean ABS      Mean dABS     Mean ddABS      Maxx ABS'')')
   10 continue
      iter=iter+1
      call blank(fil1)
      fil1(1:40)=filnam(1:40)
      fil1(ll+1:ll+18)='_SigFreeRawdat.txt'
      open(iu30,file=fil1,status='unknown')
      close(iu30,status='delete')
      open(iu30,file=fil1,status='unknown')
      call blank(fil1)
      fil1(1:40)=filnam(1:40)
      fil1(ll+1:ll+18)='_SigFreeCurves.txt'
      open(iu31,file=fil1,status='unknown')
      close(iu31,status='delete')
      open(iu31,file=fil1,status='unknown')
      call blank(fil1)
      fil1(1:40)=filnam(1:40)
      fil1(ll+1:ll+19)='_SigFreeIndices.txt'
      open(iu32,file=fil1,status='unknown')
      close(iu32,status='delete')
      open(iu32,file=fil1,status='unknown')
c
	if(rcs)then
c	   rcmax if 0 then mvf  is computed over the entire rc otherwise rcmax = length of rc to use in computing mvf
c         calculate y and store them in yrcs to compute a regional curve from them 
	   do j=1,nc
	     do i=1,mxy
	       yrcs(i,j)=-999.0
	     enddo
	     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)
	          do i=1,nreset
	             y(nss1+i)=y(nss1+i)/xm1*xbar(j)
	          enddo
	        else
	          call mnsd(y,n,xm1,sd)
	          do i=1,n
	             y(i)=y(i)/xm1*xbar(j)
	          enddo
	        endif
	     endif
	     nmax=max(n1,n)    !complete length of regional curve
	     do i=1,nmax
	       yrcs(i,j)=y(i)
	     enddo
	     iopt1(j,iter+2)=iopt
	   enddo
          call rcurve(5000,3000,nc,yrcs,xrc,f,idate,pith,flag,nmax,eq,iopt,alpha,ndt,poff,rcmax) !(xrc is the regional mean) (f is the regional curve)
	   string='Iteration       Regional Curve MVF '  
	   write(string(11:14),'(i4)')iter
	   if(plotall)call plt1(mxy,mxs,xrc,f,1,mxage,flag,string)
	   rcnt=rcnt+1
	   do jj=1,mxage
	     	     rcscrv(jj,1,rcnt)=xrc(jj);rcscrv(jj,2,rcnt)=f(jj)
	   enddo
c
	   miny=9999.
	   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
            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.5))then
                if(iopt.eq.1.or.opt.eq.2)then
                  iopt1(j,iter+2)=-66
                else
                  iopt1(j,iter+2)=6
                endif
          	  call rcurve(5000,3000,nc,yrcs,xrc,f,idate,pith,flag,nmax,eq,iopt1(j,iter+2),alpha,ndt,poff,rcmax) !xrc is the regional mean (f is the regional curve)
	   	  call plt1(mxy,mxs,xrc,f,1,mxage,flag,string)
	   	  rcnt=rcnt+1
	   	  do jj=1,mxage
	     	     rcscrv(jj,1,rcnt)=xrc(jj);rcscrv(jj,2,rcnt)=f(jj)
	   	  enddo
              endif
            endif
c
            do i=1,n
              ii=i+iskp
              wrk1(i)=y(i)
              y(i)=x(ii,j)
              if(.not.itype)then
                 sfind(ii,j)=y(i)/f(i)
                 wrk2(i)=y(i)/f(i)
              elseif(itype)then
                 wrk2(i)=y(i)-f(i)
                 sfind(ii,j)=y(i)-f(i)
                 if(wrk2(i).lt.miny)miny=wrk2(i)
              endif
            enddo
            
c optional plot current RC MVF and RCS            
	     string='Iteration #                        '  
	     write(string(13:16),'(i4)')iter
            if(plotall)call plt3(mxy,mxs,wrk2,y,f,idate,flag,nmax,j,n,ident,string,plotall)
            
            if(.not.itype)then
               call trrw(ident(j),ify,n,wrk1,nss,ifl,iu30,iuo,'mw',tit,'ssfraw',mxy,itrdb)
               call trrw(ident(j),ify,n,f,nss,ifl,iu31,iuo,'mw',tit,'ssfcur',mxy,itrdb)
               call trrw(ident(j),ify,n,wrk2,nss,ifl,iu32,iuo,'mw',tit,'ssfind',mxy,itrdb)
            endif
	   enddo 
c
          if(itype)then  !if residuals convert residuals to positive values by adding (abs(miny) + 0.1)          
             if(miny.gt.0.)then
                miny=0.000
             elseif(miny.le.0)then
                miny=abs(miny)+0.00100  !scaling factor for negative residuals
             endif
             do j=1,nc
	        ify=idate(j,1)
    	        ily=idate(j,2)
	        ily1=idate(j,3)
	        n=ily-ify+1
	        iskp=ify-ifmny
                do i=1,n
                   ii=i+iskp
                   sfind(ii,j)=sfind(ii,j)+miny
                   wrk2(i)=sfind(ii,j)
                enddo
                call trrw(ident(j),ify,n,f,nss,ifl,iu31,iuo,'mw',tit,'stdcrv',mxy,itrdb)
                call trrw(ident(j),ify,n,wrk2,nss,ifl,iu32,iuo,'mw',tit,'stdind',mxy,itrdb)
             enddo
          endif

	else  ! NO RCS
	
	   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)
                 do i=1,nreset
                    y(nss1+i)=y(nss1+i)/xm1*xbar(j)
                 enddo
              else
                 call mnsd(y,n,xm1,sd)
                 do i=1,n
                    y(i)=y(i)/xm1*xbar(j)
                 enddo
              endif
            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.5))then
                if(iopt.eq.1.or.opt.eq.2)then
                  iopt1(j,iter+2)=-66
                else
                  iopt1(j,iter+2)=5
                endif
                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,'ssfraw',mxy,itrdb)
           call trrw(ident(j),ify,n,f,nss,ifl,iu31,iuo,'mw',tit,'ssfcur',mxy,itrdb)
           call trrw(ident(j),ify,n,wrk2,nss,ifl,iu32,iuo,'mw',tit,'ssfind',mxy,itrdb)

         enddo
         
      endif   !     if RCS
c
	
c      
      write(*,'(/38x,''Mean ABS      Mean dABS     Mean ddABS      Maxx ABS'')')
c      
      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)
      call mnsd(crn2(1,1),maxn,xmstd1,sdstd1)
      do i=1,maxn
        ybar(i,1)=crn2(i,1)				!runcor bit added 01/03/11
      enddo  
      call mnsd(crn2(nss1+1,1),nreset,xm,sd)
ccc      do i=1,maxn
ccc        crn2(i,1)=crn2(i,1)-xm+1.0
ccc      enddo

      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(hpfilt.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
      testdd=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
          if(i.gt.2)then
            testdd=testdd+dabs(y(i)-2.0*y(i-1)+y(i-2))           
          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.smooth1)then
        smooth1=testd
        star2='*'
      else
        star2=' '
      endif
      testdd=testdd/dfloat(maxn1-2)
      if(testdd.lt.smooth2)then
        smooth2=testdd
        star3='*'
      else
        star3=' '
      endif
      call qcksrt(maxn1,f)
      testb=f(maxn1)
      if(testb.lt.big)then
        big=testb
        star4='*'
      else
        star4=' '
      endif
      test1(iter)=test
      test2(iter)=testd
      test3(iter)=testdd
      test4(iter)=testb
      write(*,'(1x,''Iter '',i3,''  Convergence criteria: '',f13.8,a,3(f14.8,a))')
     *  iter,test,star1,testd,star2,testdd,star3,testb,star4
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
cccccccccccccccc - New         
ccc      elseif(iter.gt.6.and.(test1(iter).gt.test1(iter-1)))then
      elseif(iter.gt.6.and.(test2(iter).gt.test2(iter-1)))then
        write(*,'(/1x,''WARNING: Convergence tests are not monotonic decreasing, iter: '',i3)')iter
        do i=1,maxn
          crn2(i,1)=crnsav(i)
        enddo
        iter=iter-1
        goto 20
cccccccccccccccc - New
      elseif(test.gt.tol)then
        do i=1,maxn
          crnsav(i)=crn2(i,1)
        enddo
        goto 10
      endif
c
   20 continue
   
      if(rcs)then
         open(unit=69,file='RCScurves.out',status='unknown')
         string='   RC     fRC '   
         write(69,'(a10,200(a15,1x))')'Bio Year ',(string,j=1,rcnt)
         do p=1,nmax
           write(69,'(400f8.2)')float(p),((rcscrv(p,j,k),j=1,2),k=1,rcnt)
         enddo
         close(69)
         call plt4(mxlen,maxit,rcscrv,nmax,rcnt)
      endif
c            
      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:23)='Original STD Chronology'
      call blank(annot2)
      annot2(1:21)='Signal Free STD Chron'
      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 *** perform running rbar variance stabilization on signal free indices
c
      if(istab.eq.'y')then
        write(*,'(/1x,''Perform signal free indices running rbar calculations'')')
        call blank(fil1)
        fil1(1:40)=filnam(1:40)
        fil1(ll+1:ll+13)='_Ssf_RBAR.txt'
        open(iu33,file=fil1,status='unknown')
        call runbar(sfind,ybar,idate,nc,rxm,runr,ncnt,nsg,nlp,mxy,mxs,iu33)
        close(iu33)
        call blank(tit)
        tit(1:40)='Signal Free Indices Running Correlations'
        call blank(xlab)
        xlab='Years'
        call blank(ylab)
        ylab='Correlations'
        call blank(annot1)
        annot1(1:14)='Running RBAR ('
        if(nsg.lt.100)then
          write(annot1(15:16),'(i2)')nsg
          annot1(17:21)=' Yrs)'
        else
          write(annot1(15:17),'(i3)')nsg
          annot1(18:22)=' Yrs)'
        endif
        call blank(annot2)
        annot2(1:20)='Constant Sample RBAR'
        ones(1)=-999.9
        call plot1(mxy,runr(1,2),runr(1,4),maxn,ifmny,xlab,ylab,tit,ones,annot1,annot2)
        do i=1,maxn
          crn3(i,2)=ybar(i,3)
        enddo
        call mnsd(crn3(1,2),maxn,xmstd2,sdstd2)
        do i=1,maxn
          datum=crn3(i,2)-xmstd2
          crn3(i,2)=datum+xmstd1
          ones=1.0
        enddo
        call blank(tit)
        tit(1:56)='Original vs. Running RBAR Stabilized Signal Free Indices'
        call blank(annot1)
        annot1(1:20)='Orig Signal Free STD'
        call blank(annot2)
        annot2(1:22)='Stable Signal Free STD'
        call plot1(mxy,crn2(1,1),crn3(1,2),maxn,ifmny,xlab,ylab,tit,ones,annot1,annot2)
      endif
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:20)='  Mean ABS Diff Test'
      call blank(annot2)
      annot2(1:19)='Mean dABS Diff Test'
      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:51)='Signal Free Chronology Iteration Differences'
      call blank(xlab)
      xlab(1:5)='Years'
      call blank(ylab)
      ylab(1:17)='Differences'
      call blank(annot1)
      annot1(1:20)='     First Iteration'
      call blank(annot2)
      annot2(1:19)='    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
ccc      close(iu29)
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr curvefits
*
      subroutine curvefits(y,f,iyr,n,n1,id,iopt,alpha,mxy,mxs)
      parameter(mxlen=5000,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
      
      do i=1,n
        if(f(i).lt.0.02)then
          f(i)=0.02
        endif
      enddo

      ier=0
      do i=1,n
        if(f(i).lt.0.02)then
          ier=1
          goto 90
        endif
      enddo
   90 if(ier.eq.1)then
        iopt=5
        call xmean(y,f,n,n1)
        return
      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
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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=5000)
      implicit real*8 (a-h,o-z)
      dimension a(mxy,4),c1(4),c2(3)
      real*4 y(n),f(n),pp
      data c1 /1.d0,-4.d0,6.d0,-2.d0/
      data c2 /0.d0,.33333333333333d0,1.33333333333333d0/
      data pct,pi / .50d0, 3.1415926535897935d0 /
      if(n .lt. 4)then
      f(1)=-9998.
      return
      endif
      if(n.gt.mxy)then
        write(*,'(/t2,''*** series too long for spline subroutine ***''/
     *           t2,'' *** enlarge mxy parameter in subroutine ***''/)')
        stop
      endif
      nm2=n-2
      v=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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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=5000)
      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
c      if(ls.eq.9)then
      if(ls.eq.8)then   !pjk
        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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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=5000)
      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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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.-999.0)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
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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=5000)
      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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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=5000)
      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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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*4 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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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)=-999.0
        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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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)=-999.0
        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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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)=-999.0
        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)=-999.0
        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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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)*88,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='mnsl'
      if(cproc.eq.'3')tag='mxsy'
      if(cproc.eq.'4')tag='comm'
      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)(85:88)='MNSL'
        if(cproc.eq.'3')itrdb(j)(85:88)='MXSY'
        if(cproc.eq.'4')itrdb(j)(85:88)='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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr arpad
*
      subroutine arpad(x,n,y,np,ls,mxy)
      real*4 x(mxy),y(mxy)
      real*8 arps(500),wk1(5000),wk2(5000),wkm(500)
c
c *** extend the series off both ends with prediction error filter
c *** forecasts and hindcasts of length ls-1 (i.e. # of lags)
c
      lstest=0.67*n
      if(lstest.lt.ls)ls=lstest
      next=ls
      np=n+2*next
      do i=1,mxy
        y(i)=0.0
      enddo
      call mnsd(x,n,xm,sd)
      do i=1,n
        ii=i+next
        y(ii)=x(i)-xm
      enddo
      iopt=0
      ip=next
      call memcof(y(next+1),n,ip,pm,arps,wk1,wk2,wkm)
      call extend(y,n,next,arps,ip,iopt)
      do i=1,np
        y(i)=y(i)+xm
      enddo
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr memcof
*
      subroutine memcof(data,n,m,pm,cof,wk1,wk2,wkm)
      implicit real*8 (a-h,o-z)
      real*4 data(1)
      real*8 cof(1),wk1(1),wk2(1),wkm(1)
      data zero/0.0d0/,one/1.0d0/,two/2.0d0/
      rn=n
      p=zero
      do 11 j=1,n
        p=p+dble(data(j))*dble(data(j))
11    continue
      pm=p/rn
      wk1(1)=dble(data(1))
      wk2(n-1)=dble(data(n))
      do 12 j=2,n-1
        wk1(j)=dble(data(j))
        wk2(j-1)=dble(data(j))
12    continue
      do 17 k=1,m
        pneum=zero
        denom=zero
        do 13 j=1,n-k
          pneum=pneum+wk1(j)*wk2(j)
          denom=denom+wk1(j)**2+wk2(j)**2
13      continue
        cof(k)=two*pneum/denom
        pm=pm*(one-cof(k)**2)
        if(k.ne.1)then
          do 14 i=1,k-1
            cof(i)=wkm(i)-cof(k)*wkm(k-i)
14        continue
        endif
        if(k.eq.m)return
        do 15 i=1,k
          wkm(i)=cof(i)
15      continue
        do 16 j=1,n-k-1
          wk1(j)=wk1(j)-wkm(k)*wk2(j)
          wk2(j)=wk2(j+1)-wkm(k)*wk1(j+1)
16      continue
17    continue
      pause 'never get here'
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr extend
*
c
c *** uses prediction error filter or autogregressive coefficients
c *** to extrapolate a time series off either or both ends into
c *** never-never-land. are you listening, peter pan or tinkerbell?
c
c *** it can be used to pad the ends of time series before applying
c *** digital filters, provide starting values for autoregressive
c *** modeling or prewhitening, and simply forecast or hindcast a
c *** time series into the past for future for whatever purpose.
c
c       y    time series of length n
c       n    number of observations
c     ncast  number of forecasts and/or hindcasts
c      phi   prediction error filter coefficients
c      ip    order of prediction error filter
c     iopt   option for forecasts and hindcasts
c              -1 hindcasts only
c               0 forecasts and hindcasts
c               1 forecasts only
c
c *** a cautionary note:
c *** this routine assumes that room is available on either
c *** end of array y to put the forecasts and/or hindcasts.
c *** thus, if hindcasts are to be made, you must have moved
c *** the observations in y foreward to make room for the ncast
c *** hindcasts. if forecasts are to be made, array y must be
c *** dimensioned long enough to hold the ncast forecasts.
c *** thus, y must be dimensioned to hold n+ncast values and,
c *** perhaps, n+2*ncast elements long in the main program,
c *** depending on the way the subroutine is used.
c
      subroutine extend(y,n,ncast,arps,ip,iopt)
      real*4 y(1)
      real*8 arps(1),pf,pb,zero
      data zero/0.0d0/
ccc      write(2,'(10f8.4)')(arps(j),j=1,ip)			!debug
      do 15 i=1,ncast
        pf=zero
        pb=zero
        nf=ncast+n+i
        nb=ncast+1-i
        do 10 j=1,ip
          if(iopt.eq.0.or.iopt.eq.+1)pf=pf+arps(j)*dble(y(nf-j))
          if(iopt.eq.0.or.iopt.eq.-1)pb=pb+arps(j)*dble(y(nb+j))
   10   continue
        y(nf)=sngl(pf)
        y(nb)=sngl(pb)
ccc        write(2,'(1x,2i5,f9.5,i4,f9.5)')i,nb,y(nb),nf,y(nf)	!debug
   15 continue
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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(5000),mny,mxy,mx,mnx	
      
      character dev*80,filename*60,ans*2
      logical bg,ok
c
      mny=100000.0;mxy=-100000.0
      red=255;green=255;blue=255
      bg=.false.
      ok=.true.;dev='Mac1'
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      
10      if(ok)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)
ccc      call pllab(xlab,ylab,tit)
      
      call pllab(xlab,ylab,' ')
      call plmtex('b',-34.5,.5,0.5,tit)
      
      
c     set line style
      call pllsty(1)
c     set line color
ccc      call plcol(15)  !color lines
      call plcol(9)  !color lines
      call plline(ny,x,wrk)            !plot wrk
ccc      call plmtex('b',-29.0,0.04,0.0,annot1)
      call plmtex('b',-32.5,0.06,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',-27.0,0.04,0.0,annot2)
      call plmtex('b',-32.5,0.55,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  
       ans='C '
	dev='Mac1'
       write(*,'(/1x,''|************* plot options menu ***************|'')')
       write(*,'(1x,''|-----------------------------------------------|'')')
       write(*,'(1x,''|       [s]ave   <[c]ontinue>   [r]eplot        |'')') 
       write(*,'(1x,''|_______________________________________________|'')')
       write(*,'(/1x,'' enter option number > '',$)')
       read(*,'(a1)')ans
       if((ans.eq.'S').or.(ans.eq.'s'))then
         dev='psc';ok=.false.;
         filename=tit
         call stringclean(filename,60)
         do p=1,60
           if(filename(p:p).eq.' ')goto 69
         enddo
69       write(filename(p:p+2),'(a3)')'.ps'
	  call plsdev(dev)
         call plsfnam(filename)
         goto 10
       elseif((ans.eq.'R').or.(ans.eq.'r'))then
         ok=.true.
         goto 10
       elseif((ans.eq.'C').or.(ans.eq.'c').or.(ans.eq.' ').or.
     *                             (ans.eq.'Y').or.(ans.eq.'y'))then
         return
	endif
c
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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=5000)
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	
      integer p
c      
      character dev*80,filename*60,ans*2
      logical bg,ok
c
      mny=100000.0;mxy=-100000.0
      red=255;green=255;blue=255
      bg=.false.;ok=.true.;dev='Mac1'
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.-999.0)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      
10      if(ok)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)
ccc      call pllab(xlab,ylab,tit)
      call pllab(xlab,ylab,' ')
      call plmtex('b',-34.5,.5,0.5,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)
ccc      call plmtex('b',-29.0,0.05,0.0,annot1)
      call plmtex('b',-32.5,0.06,0.0,annot1)

      call plcol(9)
ccc      call plmtex('b',-27.0,0.05,0.0,annot2)
      call plmtex('b',-32.5,0.55,0.0,annot2)

c
c *** end plotting
c
      call plend()   
c  
       ans='C '
	dev='Mac1'
       write(*,'(/1x,''|************* plot options menu ***************|'')')
       write(*,'(1x,''|-----------------------------------------------|'')')
       write(*,'(1x,''|       [s]ave   <[c]ontinue>   [r]eplot        |'')') 
       write(*,'(1x,''|_______________________________________________|'')')
       write(*,'(/1x,'' enter option number > '',$)')
       read(*,'(a1)')ans
       if((ans.eq.'S').or.(ans.eq.'s'))then
         dev='psc';ok=.false.;
         filename=tit
         call stringclean(filename,60)
         do p=1,60
           if(filename(p:p).eq.' ')goto 69
         enddo
69       write(filename(p:p+2),'(a3)')'.ps'
	  call plsdev(dev)
         call plsfnam(filename)
         goto 10
       elseif((ans.eq.'R').or.(ans.eq.'r'))then
         ok=.true.
         goto 10
       elseif((ans.eq.'C').or.(ans.eq.'c').or.(ans.eq.' ').or.
     *                             (ans.eq.'Y').or.(ans.eq.'y'))then
         return
	endif
c
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr runbar
*
      subroutine runbar(x,ybar,idate,nc,rxm,runr,ncnt,nsg,nlp,mxy,mxs,iu33)
c
c     x     -  input - real*4 matrix of input series being averaged
c     ybar  -  input - real*4 matrix containing a copy of the unstabilized series (*,1)
c             output - running rbar (runr) variance stabilized series, based on unsmoothed
c                     (*,2) and spline smoothed (*,3) runr series
c     idate - input - array of start and end years of input series
c     nc    - input - number of input series being averaged
c     runr  - output - real*4 array containing sample sizes, running rbars, effective sample sizes
c     nsg   -  input - running rbar window length, must be odd - incremented by one if not
c     nlp   -  input - overlap of running rbar window
c     mxy  -  input - maximum number of observations - from calling program
c     mxs  -  input - maximum number of series - from calling program
c
      parameter(maxlen=5000)
      real*4 x(mxy,mxs),ybar(mxy,3),rxm(mxy,7),runr(mxy,5),wrk1(maxlen)
      integer idate(mxs,3)
      nseg=0
      ncnt=0
      ifmn=idate(nc+1,1)
      ilmx=idate(nc+1,2)
      n=ilmx-ifmn+1   
      if(mod(nsg,2).eq.0)then
        nsg=nsg+1
        nlp=nsg-1
      endif
      nsgd2=nsg/2
      write(*,'(/1x,''Running RBAR window length: '',t37,i5/
     *           1x,''Running RBAR window overlap: '',t37,i5)')nsg,nlp
      do i=1,n,nsg-nlp
        nseg=nseg+1
        iskp=(nseg-1)*nsg-(nseg-1)*nlp
        ifskp=ifmn+iskp
        ilskp=ifskp+nsg-1
        iyr=ifskp+nsg/2
        if(ilskp.le.ilmx)then
          call rrbar(x,wrk1,idate,nc,ifskp,ilskp,nr,ave,sdev,see,nsg,mxy,mxs,maxlen)
          if(ave.ne.-99.9)then
            ncnt=ncnt+1
            if(ncnt.eq.1)then
              ndif1=iyr-idate(nc+1,1)
            endif 
            rxm(ncnt,1)=iyr
            rxm(ncnt,2)=nr
            rxm(ncnt,3)=ave
            rxm(ncnt,4)=sdev
            rxm(ncnt,5)=see            
            runr(ncnt+ndif1,2)=ave            
            rncores=0.0
            do j=1,nsg
              j1=ifskp-idate(nc+1,1)+j
              do k=1,nc
                if(x(j1,k).ne.-999.0)rncores=rncores+1.0
              enddo
            enddo
            rncores=rncores/float(nsg)
            rxm(ncnt,6)=ave/(ave+(1.0-ave)/rncores)
            rxm(ncnt,7)=rncores
          endif
        endif
      enddo      
      write(*,'(1x,''Running RBAR segments calculated: '',t37,i5)')ncnt
c
c *** pad the ends with values equal to the first and last running RBAR values
c
      do i=1,ndif1
          runr(i,2)=runr(ndif1+1,2)
      enddo
      ndif2=n-ncnt-ndif1
      do i=1,ndif2
          runr(i+ndif1+ncnt,2)=runr(ndif1+ncnt,2)
      enddo
c
c *** calculate the mean running RBARs for periods of constant sample size
c      
      icnt=1
      jcnt=0
      kcnt=0
      do i=2,n
        num1=runr(i-1,1)
        num2=runr(i,1)
        if(num2.eq.num1)then
          icnt=icnt+1
        else
          jcnt=jcnt+1
          kcnt=icnt
ccc          write(*,'(1x,''jcnt,kcnt: '',2i5)')jcnt,kcnt
          wrk1(jcnt)=kcnt
          icnt=icnt+1
        endif
        if(i.eq.n)then
          jcnt=jcnt+1
          kcnt=n
ccc          write(*,'(1x,''jcnt,kcnt: '',2i5)')jcnt,kcnt
          wrk1(jcnt)=kcnt
        endif
      enddo
      do i=1,jcnt
        sum=0.0
        if(i.eq.1)then
          kcnt1=1
        else
          kcnt1=wrk1(i-1)+1
        endif
        kcnt2=wrk1(i)
        rcnt=0.0
        do j=kcnt1,kcnt2
          rcnt=rcnt+1.0
          sum=sum+runr(j,2)
        enddo
        sum=sum/rcnt
        do j=kcnt1,kcnt2          
          runr(j,4)=sum
        enddo
      enddo      
c
c *** do the running rbar variance stabilization now
c
      call mnsd(ybar(1,1),n,xm,sd)
      do i=1,n
        runr(i,3)=runr(i,1)/(1.0+(runr(i,1)-1.0)*runr(i,2))
        runr(i,5)=runr(i,1)/(1.0+(runr(i,1)-1.0)*runr(i,4))
        iyr=ify+i-1
        datum=ybar(i,1)-xm
        ybar(i,2)=datum*sqrt(runr(i,2)*runr(i,3))+xm
        ybar(i,3)=datum*sqrt(runr(i,4)*runr(i,5))+xm
      enddo
      
      ifmn=idate(nc+1,1)
      write(iu33,'('' YEAR     NC   RBAR1   NEFF1   RBAR2   NEFF2'')')
      do i=1,n
        iyr=ifmn+i-1
        ncores=runr(i,1)
        write(iu33,'(i5,i7,4f8.3)')iyr,ncores,(runr(i,j),j=2,5)
      enddo

      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr rrbar
*
      subroutine rrbar(x,wrk1,idate,nc,ifskp,ilskp,nr,ave,sdev,see,nsg,mxy,mxs,maxlen)
      real*4 x(mxy,mxs),wrk1(maxlen)
      real*8 x1,x2,dr,rn
      integer idate(mxs,3)
      ifyf=idate(nc+1,1)
      nr=0
      x1=0.d0
      x2=0.d0
      do 20 i=1,nc-1
        ify1=idate(i,1)
        ily1=idate(i,2)
        do 10 j=i+1,nc
          ify2=idate(j,1)
          ily2=idate(j,2)
          ifmx=max0(ify1,ify2)
          ilmn=min0(ily1,ily2)
          if(ifmx.gt.ifskp)then
            goto 10
          endif
          ist=ifskp-ifyf+1
          nmin=ilmn-ifskp+1
          if(nmin.gt.nsg)nmin=nsg
          if(nmin.eq.nsg)then
            call pearsn(x(ist,i),x(ist,j),nmin,r,a,b)
            if(r.gt.0)then
              nr=nr+1
              wrk1(nr)=r
              dr=r
              x1=x1+dr
              x2=x2+dr*dr
            endif
          endif
   10   continue
   20 continue
      if(nr.lt.3)then
        ave=-99.9
        sdev=-99.9
        see=-99.9
        return
ccc      elseif(nr.eq.1)then
ccc        ave=r
ccc        sdev=0.0
ccc        see=0.0
ccc        return
      else
        rn=dble(nr)
        ave=x1/rn
        sdev=dsqrt((x2-rn*ave*ave)/(rn-1.d0))
      endif
      rnr=nr
      see=sdev/sqrt(rnr)
      rnc=nc
      pap=(rnr/(rnc*(rnc-1.0)/2.0))*100.0
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr pearsn -- pearson product-moment correlation coefficient
*
      subroutine pearsn(x,y,n,r,a,b)
      real*4 x(n),y(n)
      ax=0.
      ay=0.
      do j=1,n
        ax=ax+x(j)
        ay=ay+y(j)
      enddo
      ax=ax/n
      ay=ay/n
      sxx=0.
      syy=0.
      sxy=0.
      do j=1,n
        xt=x(j)-ax
        yt=y(j)-ay
        sxx=sxx+xt**2
        syy=syy+yt**2
        sxy=sxy+xt*yt
      enddo
      sgy=sqrt(syy/(n-1.))
      sgx=sqrt(sxx/(n-1.))
      yxr=sgy/sgx
      r=sxy/sqrt(sxx*syy)
      b=yxr*r
      a=ay-b*ax
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr intfrm
*
      subroutine intfrm(datum,frm)
      character datum*(*),frm*(*)
      l=len(datum)
      ll=0
      nb=0
      do 5 i=1,l
        if(datum(i:i).ne.' ')then
          nb=i-1
          goto 10
        endif
    5 continue
   10 do 15 i=1,l
        if(datum(i:i).ne.' ')ll=ll+1
   15 continue
      if(nb.gt.0)then
        do 20 i=1,ll
          j=i+nb
          datum(i:i)=datum(j:j)
   20   continue
        do 25 i=ll+1,l
          datum(i:i)=' '
   25   continue
      endif
      call blank(frm)
      frm(1:2)='(i'
      if(ll.lt.10)then
        write(frm(3:3),'(i1)')ll
        frm(4:4)=')'
      else
        write(frm(3:4),'(i2)')ll
        frm(5:5)=')'
      endif        
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine stringclean(string,size)
c
c remove the blank spaces in string
c
	integer size
	character string
	integer p,k
c
10	do p=1,size-1
	  if(string(p:p).eq.char(32))then
	     do k=p,size-1
	       string(k:k)=string(k+1:k+1)
	     enddo
	  endif
	enddo
	do p=1,size-1
	  if(string(p:p).eq.char(32).and.string(p+1:p+1).ne.char(32))then
	    goto 10
	  endif
	enddo
c
	return
	end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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*6
      character itrdb(3)*88
      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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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*4 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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*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)*88,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
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccc  New Subroutines for RCS Detrending  cccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine rcurve(mxy,mxs,nc,x,xrc,f,idate,pith,flag,mxage,eq,iopt,alpha,ndt,poff,rcmax)
c compute a simple age aligned population mean curve
c argument decs
	integer mxy,mxs,nc,mxage,iopt,ndt,rcmax
	real*4 alpha
	integer pith(mxy),idate(mxs,3)
	real*4 x(mxy,mxs),xrc(mxy),f(mxy),flag,eq(7)
	logical poff
c locals decs
	integer mxd, p,j,k,age,mnpoff,n1
	parameter (mxd=5000)    !must equal mxy
	real*4  num(mxd),sum(mxd)  
c	
	mxage=-9999
	do j=1,mxd
	  sum(j)=0.0;num(j)=0.0
	enddo
c
	mnpoff=1
	if(poff)mnpoff=minimum(pith,nc)
c	
	do j=1,nc
	  p=(pith(j)-mnpoff)+1
	  age=(idate(j,2)-idate(j,1))+1;
	  do k=1,age
	     sum(p)=sum(p)+x(k,j)
	     num(p)=num(p)+1.
	     p=p+1
	  enddo
	  if(age.gt.mxage)mxage=age
	enddo
c	
	do p=1,mxage
	  if(num(p).ne.0)then
	    xrc(p)=sum(p)/num(p)
	  else
	    xrc(p)=flag
	  endif
	enddo
c
	if(rcmax.eq.0)then
	    n1=mxage
	else
	    n1=rcmax
	endif
c	write(*,*)' In rcurve : rcmax= ',rcmax,' n1= ',n1,' mxage= ',mxage
c	
       call detrend2(xrc,f,mxage,n1,eq,iopt,5000,ndt,alpha)			!Hakan
c
	return
	end	
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr detrend2
* Use in conjunction with subroutine rcurve.f to perform quick and dirty regional curve
* standardization 
* Writen by ERC and only slighly modified by pjk.
*
* jdt detrending option
* alpha freidman smoother, lowess smoother
*
      subroutine detrend2(x,f,n1,n,eq,jdt,mxy,ndt,alpha)			!Hakan
c
         integer n1,n,jdt,mxy,ndt
	  real*4  x(mxy),f(mxy),eq(7)
      parameter (i1=1,i2=2,i3=3,i4=4,i5=5,i6=6,i10=10,i21=21,iu11=11)
      parameter (nemax=2500,ndmax=i5,ndmaxp1=i6,r0=0e0,r1=1e0,nmax=5000)
      real*4 pam(ndmaxp1,ndmax),yam(ndmaxp1),v(ndmax)
      real*4 span,alpha,w(nmax),sc(nmax,7)
      real*4 wk2(nmax),wk3(nmax)    !dimention must equal mxy
      real*4 rdt
c      
c
      do i=1,7
        eq(i)=0.0
      enddo
      ls=0
      ndif=0
      if(jdt.eq.-5)then
        do i=1,n
          wk2(i)=float(i)
          wk3(i)=x(i)
          w(i)=1.0
        enddo
        span=0.0
        if(n.lt.40)span=0.3
        iper=1
c        alpha=float(ndt)   !pjk
        call supsmu(n,wk2,wk3,w,iper,span,alpha,f,sc)
      elseif(jdt.eq.-4)then
        do i=1,n
          wk2(i)=x(i)
        enddo
         alpha=float(ndt)/100.  
c        rdt=float(ndt)
         call lowess(mxy,wk2,f,n,alpha)   !pjk
      elseif(jdt.eq.-3)then
        do i=1,n
          wk2(i)=x(i)
        enddo
c        do i=2,n
c          wk2(i-1)=x(i)
c        enddo
	 ndt=int(alpha)
        call msmoor(wk2,f,wk3,n,ndt,iter2,ier,mxy)
      elseif(jdt.eq.-2)then
        ls=8
        call spline1(n,x,f,ls)
      elseif(jdt.eq.-1)then
        do i=2,n
          wk2(i-1)=x(i)
        enddo
        call xmean2(x(2),n-1,xm1)
        do i=2,n
          wk3(i-1)=x(i)-x(i-1)						                !first-difference
        enddo
        x(1)=-999.0
        n=n-1
        do i=1,n
          x(i+1)=wk3(i)
        enddo
c        call xmean2(x,n,xm)
c        sum=0.0
c        do i=1,n
c          sum=sum+xm
c          f(i)=sum
c        enddo
c        call xmean2(f,n,xm2)
c        do i=1,n
c          f(i)=f(i)-xm2+xm1
c        enddo
        ls=n/3
        call spline(n,wk2,f,p,ls)
        do i=1,n
          wk2(i)=wk2(i)/f(i)
        enddo
        call moment(x(2),n,ave1,adev,sdev1,var,skew,curt,sens,ac)
        call moment(wk2,n,ave2,adev,sdev2,var,skew,curt,sens,ac)
        do i=1,n
          x(i+1)=(x(i+1)-ave1)/sdev1
          x(i+1)=x(i+1)*sdev2+ave2
        enddo
      elseif(jdt.eq.1.or.jdt.eq.2)then
        do i=1,7
          eq(i)=0.0
        enddo
c        call curve(n,x,f,eq)				!negative exponential curve
        call curve(n1,n,x,f,eq)				!Hakan
        if(eq(2).le.0.or.eq(5).lt.0)then
          do i=1,7
            eq(i)=0.0
          enddo
c          call trend(n,x,f,eq(3),eq(5))
          call trend(n1,n,x,f,eq(3),eq(5))		!Hakan
          if(jdt.eq.1)then
            jdt=4
          elseif(jdt.eq.2.and.eq(3).le.0)then
            jdt=5
          else
            do i=1,7
              eq(i)=0.0
            enddo
            jdt=5
            call xmean2(x,n,eq(5))
c            do i=1,n
            do i=1,n1
              f(i)=eq(5)					!Hakan
            enddo
          endif
        endif
      elseif(jdt.eq.3)then
        iopt=1
        do i=1,n
          wk2(i)=i
        enddo
	 nd=iopt+i2
        pam(1,1)=0.0
        pam(1,2)=0.0
        pam(1,3)=0.0
        pam(1,4)=0.0
        pam(1,5)=0.0
c        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter2,iopt,nemax,n,wk2,x,f,sumsq)
        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter2,iopt,nemax,n1,n,wk2,x,f,sumsq)		!Hakan
        slop=pam(i1,i2)
        eq(1)=pam(i1,i1)
        eq(2)=abs(pam(i1,i2))
        eq(3)=r0
        eq(5)=pam(i1,i3)
        eq(6)=r0
        if(slop.ge.0)then
          jdt=4
          do i=1,7
            eq(i)=0.0
          enddo
c          call trend(n,x,f,eq(3),eq(5))
          call trend(n1,n,x,f,eq(3),eq(5))				!Hakan
        endif
      elseif(jdt.eq.4.or.jdt.eq.5)then
c        call trend(n,x,f,eq(3),eq(5))				        !linear regression
        call trend(n1,n,x,f,eq(3),eq(5))				 	!Hakan
        if(jdt.eq.5.and.eq(3).gt.0)then
          jdt=6
          eq(3)=0.0
          eq(5)=0.0
          call xmean2(x,n,eq(5))						            !arithmetic mean
c          do i=1,n
          do i=1,n1								!Hakan
            f(i)=eq(5)
          enddo 
        endif       
      elseif(jdt.eq.6)then
        call xmean2(x,n,eq(5))						              !arithmetic mean
c        do i=1,n
        do i=1,n1								!Hakan
          f(i)=eq(5)
        enddo
      elseif(jdt.eq.7)then
        do i=1,n
          wk2(i)=i
        enddo
c        call hughdi(wk2,x,n,f,eq,flag)           !hugershoff growth curve
        call hughdi(wk2,x,n1,n,f,eq,flag)         !Hakan
        do i=1,n
          wk2(i)=i
          wk3(i)=f(i)
        enddo
        iopt=3
	       nd=iopt+i2
        pam(1,1)=eq(1)
        pam(1,2)=eq(2)
        pam(1,3)=eq(3)
        pam(1,4)=0.0
        pam(1,5)=0.0
c        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter2,iopt,nemax,n,wk2,x,f,sumsq)
        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter2,iopt,nemax,n1,n,wk2,x,f,sumsq)		!Hakan
        eq(1)=pam(i1,i1)
        eq(2)=pam(i1,i2)
        eq(3)=abs(pam(i1,i3))
        eq(5)=pam(i1,i5)
        eq(6)=pam(i1,i4)
c        if(eq(2).le.0)then
c          jdt=3
c          do i=1,7
c            eq(i)=0.0
c          enddo
c          call trend(n,x,f,eq(3),eq(5))
c        endif
        eq(7)=jdt
      elseif(jdt.eq.8)then
c        call gexp(x,n,f,wk2,wk3,eq,flag)         !general exponential curve
        call gexp(x,n1,n,f,wk2,wk3,eq,flag)       !Hakan
        do i=1,n
          wk2(i)=i
          wk3(i)=f(i)
        enddo
        iopt=2
	       nd=iopt+i2
        pam(1,1)=eq(1)
        pam(1,2)=eq(2)
        pam(1,3)=0.0
        pam(1,4)=0.0
        pam(1,5)=0.0
c        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter2,iopt,nemax,n,wk2,x,f,sumsq)
        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter2,iopt,nemax,n1,n,wk2,x,f,sumsq)		!Hakan
        eq(1)=pam(i1,i1)
        eq(2)=r1
        eq(3)=abs(pam(i1,i2))
        eq(5)=pam(i1,i4)
        eq(6)=pam(i1,i3)
        if(eq(2).le.0)then
          jdt=3
          do i=1,7
            eq(i)=0.0
          enddo
c          call trend(n,x,f,eq(3),eq(5))
          call trend(n1,n,x,f,eq(3),eq(5))			!Hakan
        endif
        eq(7)=jdt
      else
        if(jdt.gt.9)then
          ls=jdt											                      !fixed 50% frequency response
        elseif(jdt.lt.-9)then
          ls=ifix(abs((float(jdt)/100.0)*float(n)))			!%n frequency response
        endif
        do i=1,n
          wk2(i)=1.0
        enddo
ccc        call spl2w(x,n,f,wk2,p,ls,mxy)
ccc        call spl2w(mxy,n,x,wrk2,pp,f,ls)
        call spline(n,x,f,p,ls)						            !cubic spline filter
      endif
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine pithy(mxy,iu,pith,poff)
c
c argument declarations.
	integer iu,mxy,pith(mxy)
	logical poff
c local declarations
	integer p,j
	integer ifmt
	character junk*80,rfmt*40,infl*50,key*1
	logical fok
c
c begin
c
	write(*,'(//1x,''Is there a pith offset file ?'')')
	write(*,'(1x,''[t]rue-[f]alse  : '',$)')
	read(*,*)poff
	if(poff)then
10	  call blank(infl)
	  write(*,'(/1x,''Enter the name of the pith offset file.'')')
	  write(*,'(1x,''NOTE: the order of pith offset information must follow'')')
	  write(*,'(1x,''the order of the records in the raw ring width file > '')')
	  write(*,'(//1x,''Enter pith offset file name : '',$)')
	  read(*,'(a)')infl
         if(infl(1:1).eq.'/')stop
         inquire(file=infl,exist=fok)
         if(.not.fok)then
	     write(*,'(//10x,''File > '',a25,'' < not found'')')infl
	     pause'          Hit return to enter a new file name or / to quit program.'
	     goto 10
         endif
20	  write(*,'(//1x,''Select method to read the offset information.'')')
	  write(*,'(1x,''Single column offset ie. <i4> [1]'')')
	  write(*,'(1x,''User defined fortran format   [2]'')')
	  write(*,'(1x,''Enter option number  : '',$)')
	  read(*,*)ifmt
	  write(*,'(/1x,''Enter number of title/header lines'')')
	  write(*,'(1x,''to skip  [if none enter 0] : '',$,)')
	  read(*,*)p
	  select case(ifmt)
	    case(1)
	       rfmt='(i4)'
	    case(2)
		write(*,'(/1x,''Enter format to read just the '')')
		write(*,'(1x,''offset values (within parentheses)'')')
		write(*,'(1x,''[e.g.,(2x,i3)]  : '',$)')
		read(*,'(a)')rfmt
	    case default
		write(*,'(/10x,''Please enter 1 or 2 '')')
              pause'          Hit RETURN and try again.'
		goto 20
	  end select
	  open(unit=iu,file=infl,status='old')
	  write(*,'(/1x,''Reading offset file: '',a)')infl
	  do j=1,p
	     read(iu,'(a)')junk
	     write(*,'(1x,i2,''Title Line: '',a)')j,junk
	  enddo
	  poff=.true.
	  j=1	     
30	  read(iu,rfmt,end=40,err=900)pith(j)   
	     j=j+1
	     if(j.gt.mxy)goto 901
	  goto 30
	else  !no pith offset information, set all to 1
	  do j=1,mxy
	     pith(j)=1
	  end do
	  poff=.false.
	  return
	endif
40	close(iu)
c
	return
c
c error msgs
900	write(*,'(/10x,''There was an error reading the pith offset '')')
	write(*,'(10x,''values in file  : '',a)')infl
	write(*,'(10x,''Do you want to try again? or quit?'')')
	write(*,'(10x,''enter [\] to quit  : '',$)')
	read(*,*)key
	close(iu)
	if(key.ne.char(92))then
	   goto 20
	else
	   goto 1000
	endif
	goto 999
901	write(*,'(/10x,''The number of pith offset values read exceeds the'')')
	write(*,'(10x,''the dimensions of the array in which they are stored'')')

999	write(*,'(/10x,''This is a fatal error, there is no recovery '')')
1000	write(*,'(10x,''Program Terminating '')')
	stop
	end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      subroutine gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter2,iopt,nemax,ne,x,y,yfit,sumsq)
c uses the nr subroutine amoeba to find useful values for the time and
c ringwidth offsets t and k in y = a*exp(b*ln(x+t)-c*(x+t)) + k as well
c as a, b, and c.
c
c nd is the number of amoeba dimensions. it must be ² ndmax.
c pam contains the nd+1 vectors of vertex coordinants of the amoeba.
c   they are initialized below to order of magnitude and usual sign.
c yam contains the rss for each vertex.
c iter is the number of steps the amoeba needed to find the solution.
c   itmax in amoeba is currently set to ????.
c if iopt = 1, an offset gamma with b=0 will be fit.
c  "   "  " 2,  "    "     "    "   b=1   "  "   " .
c           3,  "    "     "              "  "   " .
c ne is the number of (x,y) data pairs. it must be ² nemax.
c y is the (ringwidth) data to be fit.
c x is a vector of the same length as y. it will be loaded with the
c   numbers 1 thru ne below. not the years!
c yfit is the fitted function,
c sumsq it's residual sum of squares (rss).
c
c
       subroutine gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter2,iopt,nemax,n1,ne,x,y,yfit,sumsq)	!Hakan

	parameter (i1=1,i2=2,i3=3,i4=4,i5=5,i6=6)
      
      integer ndmaxp1,ndmax,nd,iter2,iopt,nemax,ne
      
      real*4    pam(ndmaxp1,ndmax),yam(ndmaxp1),v(ndmax)
	     real*4    pami(i3,i6,i5)
	     parameter (r0=0e0,rh=-1e-2,rt=1e-1,r1=1e0)
      parameter (ftol=1e-5)
	     real*4    x(nemax),y(nemax),yfit(nemax)
c234567890         0         0         0         0         0         012
	     data pami/r0,r0,r0,r1,r1,r1,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,
     *          r0,r0,r0,r0,r0,r0,rh,rh,r1,r0,r0,r0,r0,r0,r0,r0,r0,r0,
     *          r0,r0,r0,r0,r0,r0,r0,r0,r0,rt,r1,rh,r0,r0,r0,r0,r0,r0,
     *          r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,rt,r1,r0,r0,r0,
     *          r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,rt/
c
c initialize pam and yam.
c the first vertex has been passed in pam as input.
c replace the original(0) vertex with this linearized solution.
      do ic=i1,i5
	       v(ic)=pam(i1,ic)
	     enddo
	     yam(i1)=gamoff(v,ndmax,nemax,ne,x,y,iopt)
      do ir=i2,nd+i1
	       do ic=i1,nd
		        pamu=pami(iopt,ir,ic)
		        pam(ir,ic)=pamu
		        v(ic)=pamu
		      enddo
	       yam(ir)=gamoff(v,ndmax,nemax,ne,x,y,iopt)
	     enddo
c
c call amoeba.
    2 call amoeba(pam,yam,ndmaxp1,ndmax,nd,iter2,ftol,
     *            iopt,nemax,ne,x,y)
	     sumsq=yam(i1)
c
c get the fitted function.
	     go to (21,22,23) iopt
   21 a=pam(i1,i1)
	     b=r0
      c=pam(i1,i2)
      t=r0
      rk=pam(i1,i3)
	     go to 9
   22 a=pam(i1,i1)
	     b=r1
	     c=pam(i1,i2)
	     t=pam(i1,i3)
	     rk=pam(i1,i4)
	     go to 9
   23 a=pam(i1,i1)
   	  b=pam(i1,i2)
   	  c=pam(i1,i3)
   	  t=pam(i1,i4)
   	  rk=pam(i1,i5)
c    9 do ie=i1,ne
    9 do ie=i1,n1						!Hakan
	       time=float(ie)+t
	       yfit(ie)=a*exp(b*alog(time)+c*time)+rk
	     enddo
	     return
	     end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	function gamoff(v,ndmax,nemax,ne,x,y,iopt)
c will compute y = a*exp(b*ln(x+t)-c*(x+t)) + k for x=1,ne.
c
      parameter (i1=1,i2=2,i3=3,i4=4,i5=5)
	     real*4    v(ndmax),x(nemax),y(nemax)
	     real*8    ad,bd,cd,td,rkd,sumsq,xd,tmd,d0,d1,yd,difd
	     parameter (d0=0d0,d1=1d0)
c
c dimensions...
      if(ne.gt.nemax) then
	       write(*,'(/1x,''in gamoff ne > nemax: '',2i10)') ne,nemax
		      pause'bye.'
		      stop
	     endif
c
c compute the residuals.
	     if(iopt.eq.i1) then
	       ad=dble(v(i1))
	       bd=d0
		      cd=dble(v(i2))
		      td=d0
		      rkd=dble(v(i3))
	       go to 8
	     endif
	     if(iopt.eq.i2) then
	       ad=dble(v(i1))
	       bd=d1
		      cd=dble(v(i2))
		      td=dble(v(i3))
		      rkd=dble(v(i4))
	       go to 8
	     endif
	     if(iopt.eq.i3) then
	       ad=dble(v(i1))
	       bd=dble(v(i2))
		      cd=dble(v(i3))
		      td=dble(v(i4))
		      rkd=dble(v(i5))
	       go to 8
	     endif
    8 sumsq=d0
	     do ie=i1,ne
	       xd=dble(x(ie))
		      tmd=xd+td
		      if(tmd.lt.d0) then
		        write(*,'(/1x,''in gamoff tmd ² 0 at xd.'',2d17.10)') tmd,xd
c23456789012345678901234567890123456789012345678901234567890123456789012
          pause'this is a fatal condition.'
		        stop
		      endif
		      yd=ad*dexp(bd*dlog(tmd)+cd*tmd)+rkd
		      difd=yd-dble(y(ie))
		      sumsq=sumsq+difd*difd
	     enddo
	     gamoff=sngl(sumsq)
	     return
	     end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine amoeba(p,y,mp,np,ndim,iter2,ftol,iopt,
     *                 nemax,ne,x,ye)
      integer iter2,mp,ndim,np,nmax,itmax
      real*4 ftol,p(mp,np),y(mp)
	     real*4    x(nemax),ye(nemax)
      parameter (nmax=5,itmax=2000)
cu    uses ametry,gamoff
      integer i,ihi,ilo,inhi,j,m,n
      real*4 rtol,sum,swap,ysave,ytry,psum(nmax)
      iter2=0
1     do 12 n=1,ndim
        sum=0.
        do 11 m=1,ndim+1
          sum=sum+p(m,n)
11      continue
        psum(n)=sum
12    continue
2     ilo=1
      if (y(1).gt.y(2)) then
        ihi=1
        inhi=2
      else
        ihi=2
        inhi=1
      endif
      do 13 i=1,ndim+1
        if(y(i).le.y(ilo)) ilo=i
        if(y(i).gt.y(ihi)) then
          inhi=ihi
          ihi=i
        else if(y(i).gt.y(inhi)) then
          if(i.ne.ihi) inhi=i
        endif
13    continue
      rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo)))
      if (rtol.lt.ftol) then
        swap=y(1)
        y(1)=y(ilo)
        y(ilo)=swap
        do 14 n=1,ndim
          swap=p(1,n)
          p(1,n)=p(ilo,n)
          p(ilo,n)=swap
14      continue
        return
      endif
      if (iter2.ge.itmax) then
        pause 'itmax exceeded in amoeba'
        return
      endif
      iter2=iter2+2
c	  write(*,'(1x,''iter2,y: '',i5,3e14.7)') iter2,(y(iy),iy=1,3)
c	  write(*,'(1x,''p: '')')
c	  do iv=1,3
c	    write(*,'(1x,2e14.7)') (p(iv,itk),itk=1,2)
c	  enddo
c	  write(*,'(1x,''a,b,c: '',3e14.7)') cc(1),cc(2),cc(3)
c234567890         0         0         0         0         0         012
      ytry=ametry(p,y,psum,mp,np,ndim,ihi,-1.0,iopt,
     *            nemax,ne,x,ye)
      if (ytry.le.y(ilo)) then
        ytry=ametry(p,y,psum,mp,np,ndim,ihi,2.0,iopt,
     *              nemax,ne,x,ye)
      else if (ytry.ge.y(inhi)) then
        ysave=y(ihi)
        ytry=ametry(p,y,psum,mp,np,ndim,ihi,0.5,iopt,
     *              nemax,ne,x,ye)
        if (ytry.ge.ysave) then
          do 16 i=1,ndim+1
            if(i.ne.ilo)then
              do 15 j=1,ndim
                psum(j)=0.5*(p(i,j)+p(ilo,j))
                p(i,j)=psum(j)
15            continue
			  itag=8
c234567890         0         0         0         0         0         012
              y(i)=gamoff(psum,np,nemax,ne,x,ye,iopt)
            endif
16        continue
          iter2=iter2+ndim
          goto 1
        endif
      else
        iter2=iter2-1
      endif
      goto 2
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      function ametry(p,y,psum,mp,np,ndim,ihi,fac,iopt,
     *                nemax,ne,xe,ye)
      integer ihi,mp,ndim,np,nmax
      real*4 fac,p(mp,np),psum(np),y(mp)
	     real*4    xe(nemax),ye(nemax)
      parameter (nmax=5)
cu    uses gamoff
      integer j
      real*4 fac1,fac2,ytry,ptry(nmax)
      fac1=(1.-fac)/ndim
      fac2=fac1-fac
      do 11 j=1,ndim
        ptry(j)=psum(j)*fac1-p(ihi,j)*fac2
11    continue
      i1=1
	     ndimm1=ndim-i1
	     r0=0e0
      if(iopt.gt.i1) then
		      if(ptry(ndimm1).le.r0) ptry(ndimm1)=r0
	     endif
	     itag=9
      ytry = gamoff(ptry,np,nemax,ne,xe,ye,iopt)
      if (ytry.lt.y(ihi)) then
        y(ihi)=ytry
        do 12 j=1,ndim
          psum(j)=psum(j)-p(ihi,j)+ptry(j)
          p(ihi,j)=ptry(j)
12      continue
      endif
      ametry=ytry
      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr hughdi
*
c      subroutine hughdi(x,y,ntot,xy,coef,flag)
      subroutine hughdi(x,y,n1,ntot,xy,coef,flag)		!Hakan
c
c     x     --  predictor variable -- time t=1,ntot
c     y     --  tree-ring series
c     ntot  --  number of years
c     xy    --  hugershoff growth curve estimates
c     coef  --  vector of hugershoff coefficients
c     flag  --  missing value flag
c
c     ausgleich mit hugershoff-funktion
c     schaetzen der parameter mit direkter berechnung (h.p. holzhauser)
c     a, cm, ck    im modell 
c     y(i) = a * t(i)**cm * exp(ck*t(i))   fuer i=1....ntot
c     fuer die berechnung muessen folgender summenregister vom 
c     logarithmierten x und y und unlogarithmierten x (t) vorbereitet werden.
c     sumt, sumt2, sumx, sumx2, sumy, sumyt, sumxy und sumxt.
c     regression ueber hilfsfunktion 
c     z = ln(y) = ln(a) + m*ln(t) - k*t 
c     umgesetzt in z = a' + m*t' - k*t
c     nrow entsprechend dimensionierung hauptprogramm
c     ncoef muss 3 sein
c
c     september 1989 / nos/ve /  braeker/nogler
c     pc-version 2.1 / juni 1992 / p. nogler
c     modified by e.r. cook oct 23, 1996
c
      implicit real*8 (a-h,o-z)
      real*4 x(ntot),y(ntot),xy(ntot),coef(7),flag
      real*8 a,ck,cm,se,sf,cc,dcoef(3)
      integer ntot,ncoef,nred
      logical lber
      lber=.true.
      cc = -999.0d0
c
c     beginn hughdi
c
      if(.not.lber) goto 500 
      nred  = 0 
      sumt  = 0.0d0
      sumx  = 0.0d0
      sumx2 = 0.0d0
      sumt2 = 0.0d0
      sumxt = 0.0d0
      sumy  = 0.0d0
      sumxy = 0.0d0
      sumyt = 0.0d0
c
      do 200 k = 1,ntot
        xk=dble(x(k))
        yk=dble(y(k))
        if(yk.eq.flag) goto 200
        if(yk.le.0.0d0) goto 200
        if(xk.le.0.0d0) goto 200
        nred = nred+1 
        xl = dlog(xk)
        yl = dlog(yk)
        sumt  = sumt  + xk
        sumt2 = sumt2 + xk*xk
        sumx  = sumx  + xl
        sumx2 = sumx2 + xl*xl
        sumxt = sumxt + xl*xk
        sumy  = sumy  + yl
        sumxy = sumxy + xl*yl
        sumyt = sumyt + yl*xk
  200 continue
c 
      rn1 = (sumx2 * sumt2 - sumxt**2) 
      rn2 = (sumt * sumxt - sumx * sumt2) 
      rn3 = (sumx * sumxt - sumt * sumx2)
      rn  = nred * rn1 + sumx * rn2 + sumt * rn3
      rza = sumy * rn1 + sumxy * rn2 + sumyt * rn3
      a   = rza / rn
      rz2 = (sumx * sumt - nred * sumxt)
      rzm = sumy * rn2 + sumxy * (nred * sumt2 - sumt**2) + sumyt * rz2
      cm  = rzm / rn
      rz1 = (sumt * sumx2 - sumx * sumxt)
      rz2 = (nred * sumxt - sumx * sumt)
      rzk = sumy * rz1 + sumxy * rz2 + sumyt * (sumx**2 - nred * sumx2)
      ck  = rzk / rn
c
      ncoef = 3
      dcoef(1) = dexp(a)
      dcoef(2) = cm
      dcoef(3) = -ck
      coef(1) = sngl(dexp(a))
      coef(2) = sngl(cm)  
      coef(3) = sngl(-ck) 
ccc      if(coef(2).le.0.0)then
ccc        coef(2)=-999.0
ccc        xy(1)=-999.0
ccc        goto 700
ccc      endif
c 
c     coef-korrekturfaktor berechnen ueber abweichungen
c
      se = 0.0d0
      sf = 0.0d0
c      do k = 1,ntot 
      do k = 1,n1							!Hakan 
        xk = dble(x(k))
        yk = dble(y(k))
        dxy = dcoef(1) * xk**dcoef(2) * dexp(dcoef(3)*xk)
        xy(k) = sngl(dxy)
        if(yk.ne.flag)then 
          se = se+(yk-dxy) 
          sf = sf+dxy 
        endif
      enddo
      cc = se/sf+1.0d0
      if(cc.eq.0.0d0) goto 700 
      dcoef(1) = dcoef(1)*cc
      coef(1) = sngl(dcoef(1))
c      do k = 1,ntot 
      do k = 1,n1							!Hakan 
        xy(k) = xy(k)*cc 
      enddo
      goto 700
  500 continue
c     berechnen der ausgleichskurve und xy in kolonne 1 setzen
c      do k = 1,ntot 
      do k = 1,n1							!Hakan
        xk = dble(x(k))
        dxy = dcoef(1) * xk**dcoef(2) * dexp(dcoef(3)*xk)
        xy(k) = sngl(dxy) 
      enddo
  700 return
      end 
*sbr gexp
*
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine gexp(y,nt,n,f,wk1,wk2,coef,flag)		!Hakan
c
c     y     --  tree-ring series
c     n     --  number of years
c     f     --  general exponential curve estimates
c     wk1   --  work vector
c     wk2   --  work vector
c     coef  --  vector of hugershoff coefficients
c     flag  --  missing value flag
c
      real*4 y(n),f(n),wk1(n),wk2(n),coef(7),flag
      real*8 a,b,r,ax,ay,dn1,xt,yt,sxx,syy,sxy,sgy,sgx,yxr,ri
      if(n.le.0)then
        write(*,'(/1x,''*** n .le. 0 in sbr gexp ***''/)')
        pause
        stop
      endif
      a=0.0d0
      b=0.0d0
      n1=0
c
c *** linearize the problem
c
      do i=1,n
        if(y(i).gt.0.and.y(i).ne.flag)then
          n1=n1+1
          wk1(n1)=float(i)
          wk2(n1)=alog(y(i))-alog(wk1(i))
        endif
      enddo
      if(n1.le.0)then
        write(*,'(/1x,''*** n1 .le. 0 in sbr gexp ***''/)')
        pause
        stop
      endif
c
c *** now compute simple linear regression
c
      dn1=dfloat(n1)
      ax=0.0d0
      ay=0.0d0
      do j=1,n1
        ax=ax+dble(wk1(j))
        ay=ay+dble(wk2(j))
      enddo
      ax=ax/dn1
      ay=ay/dn1
      sxx=0.0d0
      syy=0.0d0
      sxy=0.0d0
      do j=1,n1
        xt=dble(wk1(j))-ax
        yt=dble(wk2(j))-ay
        sxx=sxx+xt**2
        syy=syy+yt**2
        sxy=sxy+xt*yt
      enddo
      sgy=dsqrt(syy/(dn1-1.0d0))
      sgx=dsqrt(sxx/(dn1-1.0d0))
      yxr=sgy/sgx
      r=sxy/dsqrt(sxx*syy)
      b=yxr*r
      a=ay-b*ax
c
c *** unlinearize the fitted model and calculate growth curve values
c
      a=dexp(a)
c      do i=1,n
      do i=1,nt					!Hakan
        ri=dble(i)
        f(i)=sngl(a*ri*dexp(b*ri))
      enddo
      coef(1)=sngl(a)
      coef(2)=-sngl(b)
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr moment
*
      subroutine moment(data,n,ave,adev,sdev,var,skew,curt,sens,ac)
      real*4 data(n)
      if(n.eq.0)then
        ave=-99.0
        adev=-99.0
        sdev=-99.0
        var=-99.0
        skew=-99.0
        curt=-99.0
        sens=-99.0
        ac=-99.0
        return
      endif
      rn=n
      s=0.0
      do j=1,n
        s=s+data(j)
      enddo
      ave=s/rn
      if(n.eq.1)then
        adev=-99.0
        sdev=-99.0
        var=-99.0
        skew=-99.0
        curt=-99.0
        sens=-99.0
        ac=-99.0
        return
      endif
      adev=0.0
      var=0.0
      skew=0.0
      curt=0.0
      sens=0.0
      do j=1,n
        s=data(j)-ave
        adev=adev+abs(s)
        p=s*s
        var=var+p
        p=p*s
        skew=skew+p
        p=p*s
        curt=curt+p
        if(j.lt.n)then
          rnum=data(j+1)-data(j)
          dnum=(data(j+1)+data(j))/2.0
          if(dnum.gt.0)then
            sens=sens+abs(rnum/dnum)
          endif
        endif
      enddo
      sens=sens/(rn-1.0)
      adev=adev/rn
      var=var/(rn-1.0)
      bvar=var*(rn-1.0)/rn
      if(var.gt.0.0)then
        sdev=sqrt(var)
        bsdev=sqrt(bvar)
        skew=skew/(rn*sdev**3)
        curt=curt/(rn*var**2)
        if(n.lt.3)then
          bfs=1.0
        else
          bfs=(sqrt(abs((rn*(rn-1.0))))/(rn-2.0))*(1.0+8.5/rn)
        endif
        if(n.lt.4)then
          bfk=1.0
        else
          bfk=(rn*rn*(rn+1.0))/((rn-1.0)*(rn-2.0)*(rn-3.0))
        endif
        skew=skew*bfs
        curt=curt*bfk
        do i=1,n
	         data(i)=(data(i)-ave)/bsdev
        enddo
        cxx=0.0
        do i=2,n
          cxx=cxx+data(i)*data(i-1)
        enddo
        ac=cxx/rn
        do i=1,n
          data(i)=data(i)*bsdev+ave
        enddo
      else
        sdev=-99.0
        skew=-99.0
        curt=-99.0
        ac=-99.0
      endif
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr xmean2
*
      subroutine xmean2(x,n,xm)
      real*4 x(n)
      if(n.le.0)then
        write(*,'(/1x,''error in xmean2: n=0'')')
        pause
        stop
      endif
      rn=n
      xm=0.0
      do i=1,n
        xm=xm+x(i)
      enddo
      xm=xm/rn
      return
      end

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

         goto  1
c rightmost pt (may be greater than nright because of ties)
   7  nrt = j-1
      if (a .gt. 0.0) goto 8
         ok = .false.
         goto  16
   8     ok = .true.
c weighted least squares
         do  9 j = nleft, nrt
c make sum of w(j) == 1
            w(j) = w(j)/a
   9        continue
         if (h .le. 0.) goto 14
            a = 0.0
c use linear fit
            do  10 j = nleft, nrt
c weighted center of x values
               a = a+w(j)*x(j)
  10           continue
            b = xs-a
            c = 0.0
            do  11 j = nleft, nrt
               c = c+w(j)*(x(j)-a)**2
  11           continue
            if (sqrt(c) .le. .001*range) goto 13
               b = b/c
c points are spread out enough to compute slope
               do  12 j = nleft, nrt
                  w(j) = w(j)*(b*(x(j)-a)+1.0)
  12              continue
  13        continue
  14     ys = 0.0
         do  15 j = nleft, nrt
            ys = ys+w(j)*y(j)
  15        continue
  16  return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr msmoor
*
      subroutine msmoor(x,y,z,n,m,iter2,ier,mxy)
c     repeatedly calls msmoo until either absum < small or iter2 > itmax.
c
      parameter (i0=0,i1=1,i2=2,i3=3)
      parameter (itmax=10,small=1e-6)                                                                               
      real*4  x(mxy),y(mxy),z(mxy)
c
c input error checks
      ier=i0
      if(n.lt.m.or.n.gt.mxy)ier=i1
      if(m.lt.i1.or.m.gt.mxy)ier=i2
      if(mod(m,i2).ne.i1)m=m-1
      if(mod(m,i2).ne.i1)ier=i3
      if(ier.ne.i0) then
        write(*,'(/1x,''input error '',i1,'' in msmoor.'')') ier
        pause
        return
      endif
c pass.
      iter2=i0
   10 iter2=iter2+i1
      if(iter2.gt.itmax) then
        ier=20
		      write(*,'(/1x,''*** maximum iterations reached in '',
     *    ''msmoor ***'')')
        return
      endif
      call msmoo(x,y,z,n,m,iter2,ier,mxy)
      if(ier.ne.i0) then
        ier=ier+10
        return
      endif
      absum=0e0
      do i=i1,n                                                              
        absum = absum + abs(x(i) - y(i))
      enddo
c      write(*,'(/1x,i5,f10.7)') iter2,absum
c      write(11,'(/1x,i5,f10.7)') iter2,absum
      if(absum.lt.small) return
      do i=i1,n
        x(i)=y(i)
      enddo
      go to 10
      end 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr msmoo
*                                                                      
      subroutine msmoo(x,y,z,n,m,iter2,ier,mxy)
c     will do median smoothing of a time series. n is the number of
c     data in the series, i.e. the series length. m is the filter
c     length which must be odd. x comes in,
c                               y goes out, and
c                               z is a little work vector.
c                               oh la di da!
c     at the ends medians of 1, 3, 5,..., m-2 data points are taken.
c
      parameter (i0=0,i1=1,i2=2)
      parameter (r0=0e0)
      real*4 x(mxy),y(mxy),z(mxy)
c
c
c input error checks
      ier=i0
      if(n.lt.m.or.n.gt.mxy) ier=i1
      if(m.lt.i1.or.m.gt.mxy) ier=i2
      if(mod(m,i2).ne.i1) ier=3
      if(ier.ne.i0) then
        write(*,'(/1x,''input error '',i1,'' in msmoo.'')') ier
        pause
        return
      endif
      if(m.eq.i1) then
        do i=i1,n
          y(i)=x(i)
        enddo
        return
      endif
      mm = (m-i1)/i2
c
c save the first point.
ccc      y(i1)=x(i1+1)
      y(i1)=x(i1)
c
c do from 2 thru mm.
      if(mm.gt.i1) then
        do i=i2,mm
          im=i2*i-i1
          do j=i1,im
            z(j) = x(j)
          enddo
ccc          call rank(z,im,mmax)
          call sort(im,z)
          y(i) = z(i)
        enddo
      endif
c
c do the middle.
      do i=mm+i1,n-mm
        jj=i-mm-i1
        do j=i-mm,i+mm
          z(j-jj) = x(j)
        enddo
ccc        call rank(z,m,mmax)
        call sort(m,z)
        y(i) = z(mm+i1)
      enddo
c
c do from n-mm+1 thru n-1.
      if(mm.gt.i1) then
        np1=n+i1
        do i=i2,mm
          im=i2*i-i1
          do j=i1,im
            z(j) = x(np1-j)
          enddo
ccc          call rank(z,im,mxy)
          call sort(im,z)
          y(np1-i)=z(i)
        enddo
      endif
c
c save the last point.
ccc      y(n)=x(n-1)
      y(n)=x(n)
c
      do i=i1,mm
        j=n-i+1
        y(i)=y(mm+1)
        y(j)=y(n-mm)
      enddo
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr sort
*
*      heapsort algorithm from numerical recipes
*
*      n   =  number of observations in ra to be sorted
*      ra  =  input series to be sorted. ra is returned
*             with its observations ranked in ascending
*             numerical order.
*
      subroutine sort(n,ra)
      real*4 ra(n)
      l=n/2+1
      ir=n
   10 continue
        if(l.gt.1)then
          l=l-1
          rra=ra(l)
        else
          rra=ra(ir)
          ra(ir)=ra(1)
          ir=ir-1
          if(ir.eq.1)then
            ra(1)=rra
            return
          endif
        endif
        i=l
        j=l+l
   20   if(j.le.ir)then
          if(j.lt.ir)then
            if(ra(j).lt.ra(j+1))j=j+1
          endif
          if(rra.lt.ra(j))then
            ra(i)=ra(j)
            i=j
            j=j+j
          else
            j=ir+1
          endif
        go to 20
        endif
        ra(i)=rra
      go to 10
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c function to find the maximum value
	function minimum(iarray, nsize)
	dimension iarray(nsize)
	minimum=iarray(1)
	do 10 j=2,nsize
	   if (minimum.gt.iarray(j)) minimum=iarray(j)
10     continue
	return
	end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c function to find the maximum value
	function maximum(iarray, nsize)
	dimension iarray(nsize)
	maximum=iarray(1)
	do 10 j=2,nsize
	   if (maximum.lt.iarray(j)) maximum=iarray(j)
10     continue
	return
	end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr 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)
      integer mxs,ns
      integer n(mxs),iyr(mxs)
      character*8 id(mxs)
c
c local decs
	integer mxs2,i,l,m,j,k,ifl,icy,lcy,nly
	integer lyr,iz,lly,is,ici,nyc,ny,nci,nsci,nrci
	parameter (mxs2=3000)   !must = mxrec
       integer l20(mxs2,20)
	character sqc(mxs2)*40,sa*40,sb*40
c
*
      if(ns.gt.mxs2)goto 100
      ifl=0
c
      icy= 32000
      lcy=-32000
*
* here is where you would dump your 
* file data into the character array sqc 
      do 10 m=1,ns
          lyr=iyr(m)+n(m)-1
          icy=min(iyr(m),icy)
          lcy=max(lyr,lcy)
          iz=9000-iyr(m)
c          write(sqc(ns),'(5i5,2x,a8,5x)')lyr,iz,n(m),iyr(m),ns,id(m)
          write(sqc(m),'(5i5,2x,a8,5x)')lyr,iz,n(m),iyr(m),ns,id(m)
c          write(*,'(5i5,2x,a8,5x)')lyr,iz,n(m),iyr(m),ns,id(m)
c          pause '5'
10 	continue
c
       do 20 i=1,ns
          sa=sqc(i)
          sb=sqc(i)
          k=i
          do 30 j=i,ns
             if(sqc(j) .gt. sa)then
                sa=sqc(j)
                k=j
             endif
   30     continue
          sqc(k)=sb
20     sqc(i)=sa
*
* make array of starting and ending years of series
       nly=1
c            pause'1-commint'
       read(sqc(1),'(i5,t1,i5,t16,i5)')l20(1,3),l20(1,2),l20(1,1)
       lly=l20(nly,3)
       do 60 i=2,ns
          read(sqc(i),'(i5,t1,i5,t16,i5)',err=1000)l20(mxs,3),l20(i,2),l20(i,1)
          if(lly .ne. l20(mxs,3))then
             nly=nly+1
             l20(nly,3)=l20(mxs,3)
             lly=l20(nly,3)
          endif
60     continue
c	write(*,*)nly
c	pause '2-commint'
       nrci=0
       do 40 l=1,nly
          do 50 i=1,ns
             nc=0
             ny=l20(l,3)-l20(i,1)+1
             do 70 is=1,ns
                 if(l20(is,1) .le. l20(i,1) .and. l20(is,2) .ge. l20(l,3))then
                    nc=nc+1
                 endif
70           continue
             nyc=nc*ny
             if(nyc .gt. nrci)then
                 nrci=nyc
                 ici=l20(i,1)
                 lci=l20(l,3)
                 nci=ny
                 nsci=nc
             endif
50        continue
c	   write(*,*)l,nly,i,ns
c	   pause '3-commint'
40     continue
c	pause '4-commint'
*
* print results on screen
       write(*,'(/3x,42(''*''))')
       write(*,'(2x,'' **  Time span of chronology: '',i5,'' to'',i5)')
     +        icy,lcy
       write(*,'(2x,'' **  length =               : '',i13)')
     +        lcy-icy+1
       write(*,'(2x,'' **  Number of series in chronology: '',i6)')ns
       write(*,'(2x,'' **  Optimum common interval: '',i5,'' to'',i5)')
     +        ici,lci
       write(*,'(2x,'' **  Number of years included: '',i12)')nci
       write(*,'(2x,'' **  Number of Series included:     '',i7,/,
     + 2x,'' **  Number of rings:               '',i7)')nsci,nrci
       write(*,'(2x,'' **  Length of longest Series: '',i12)')maximum(n,ns)
       write(*,'(2x,'' **  Length of shortest Series: '',i11)')minimum(n,ns)

       write(*,'(3x,42(''*''))')
       return
100    write(*,'('' >> file contains more than program limit of'',i4,
     +'' series'')')mxs2-1
       stop' too many series for subroutine Commint'
1000   write(*,*)sqc(i)
	write(*,'(1x,''error reading sqc)'')')
	return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine plt1(mxsl,rmax,mey,ycal,fycaly,lycaly,flag,string)
c
c	argument
	integer mxsl,rmax
	integer fycaly,lycaly
	real*4 mey(mxsl),ycal(mxsl)
	character string*30
c	local
	integer k
	integer red,blue,green
	real*4 miny, maxy
	real*4 wrk1(5000)
	character ans*1,filename*9,dev*80
	logical ok,bg
c
       dev='Mac1';ok=.true.
       red=255;green=255;blue=255
       bg=.false.
c
10	if(ok)call plsdev(dev)
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	
	miny=9999.;maxy=-9999.
	do k=fycaly,lycaly
	    if(mey(k).ne.flag.and.mey(k).gt.maxy)maxy=mey(k)
	    if(mey(k).ne.flag.and.mey(k).lt.miny)miny=mey(k)
	    if(ycal(k).ne.flag.and.ycal(k).lt.miny)miny=ycal(k)
	    if(ycal(k).ne.flag.and.ycal(k).gt.maxy)maxy=ycal(k)
	    wrk1(k)=real(k)
	enddo
	if(miny.ge.0.)miny=0.
c
c initialize plplot
	call plinit()
c set lable color
	call plcol(15)   !color box and lables
c define the viewport one column by one rows
	call plssub(1,1)
	call plenv(real(fycaly),real(lycaly),miny,maxy,0,2)
	call pllab('Biological Year','Y',string)
c	step=-12.
	call plcol(1)   !color box and lables
c	write(tit(1:20),'(i5)')nboot
c	call plmtex('rv',step,1.07,0.0,tit)
c	step=-6.
c	call plmtex('rv',step,1.07,0.0,'Reps')
	call plline(lycaly,wrk1,mey)
	call plcol(9)   !color box and lables
	call plwid(2)
	call plline(lycaly,wrk1,ycal)
	call plend()
c	
       ans='C '
	dev='Mac1'
       write(*,'(/1x,''|************* plot options menu ***************|'')')
       write(*,'(1x,''|-----------------------------------------------|'')')
       write(*,'(1x,''|       [s]ave   <[c]ontinue>   [r]eplot        |'')') 
       write(*,'(1x,''|_______________________________________________|'')')
       write(*,'(/1x,'' enter option number > '',$)')
       read(*,'(a1)')ans
       if((ans.eq.'S').or.(ans.eq.'s'))then
         dev='psc';ok=.false.;filename='fig11.ps'
	  call plsdev(dev)
         call plsfnam(filename)
         goto 10
       elseif((ans.eq.'R').or.(ans.eq.'r'))then
         ok=.true.
         goto 10
       elseif((ans.eq.'C').or.(ans.eq.'c').or.(ans.eq.' ').or.
     *                             (ans.eq.'Y').or.(ans.eq.'y'))then
         return
	endif
c
	return
	end
c

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine plt2(mxy,mxs,pltdata,nc,idate,flag,string)
c
c	argument
	integer mxy,mxs,nc
	integer idate(mxs,3)
	real*4 pltdata(mxy,mxs),flag
	character string*30
c	local
	integer p,j,k, minx, maxx
	integer red,blue,green
	real*4 miny, maxy
	real*4 wrk1(5000),wrk2(5000),wrk3(5000),wrk4(5000)
	character ans*1,filename*9,dev*80
	logical ok,bg
c
       dev='Mac1';ok=.true.
       red=255;green=255;blue=255
       bg=.false.
c
10	if(ok)call plsdev(dev)
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	
	miny=9999.;maxy=-9999.;minx=9999;maxx=-9999
	do p=1,nc
	  j=1
	  do k=idate(p,1),idate(p,2)
	      if(pltdata(j,p).ne.flag.and.pltdata(j,p).gt.maxy)maxy=pltdata(j,p)
	      if(pltdata(j,p).ne.flag.and.pltdata(j,p).lt.miny)miny=pltdata(j,p)
	      j=j+1
	  enddo
	  if(idate(p,1).lt.minx)minx=idate(p,1)
	  if(idate(p,2).gt.maxx)maxx=idate(p,2)
	enddo
c
c initialize plplot
	call plinit()
c set lable color
	call plcol(15)   !color box and lables
c define the viewport one column by one rows
	call plssub(1,2)
	call plenv(real(minx),real(maxx),miny,maxy,0,2)
	call pllab('Calendar Year','Y','First Detrending')
	call plcol(1)   !color box and lables
	do p=1,nc
	   k=1
	   do j=idate(p,1),idate(p,2)
	     wrk1(k)=real(j);wrk2(k)=pltdata(k,p)
	     k=k+1
	   enddo
	   k=k-1
	    call plline(k,wrk1,wrk2)
	enddo
c
c set lable color
	call plcol(15)   !color box and lables
	call plenv(real(minx),real(maxx),miny,maxy,0,2)
	call pllab('Calendar Year','Y','Average Indices')
	do p=1,mxy
	   wrk1(p)=0.;wrk2(p)=0.;wrk3(p)=0
	enddo
	do p=1,nc
	   k=1
	   do j=idate(p,1),idate(p,2)
	     wrk2(j)=wrk2(j)+pltdata(k,p)
	     wrk3(j)=wrk3(j)+1.
	     k=k+1
	   enddo
	   k=k-1
	enddo
	k=1
	do j=minx,maxx
	   if(wrk3(j).ne.0)wrk1(k)=wrk2(j)/wrk3(j)
	   wrk4(k)=real(j)
	   k=k+1
	enddo
      	k=k-1
	call plcol(1)   !color box and lables
	call plline(k,wrk4,wrk1)
	call plend()
c	
       ans='C '
	dev='Mac1'
       write(*,'(/1x,''|************* plot options menu ***************|'')')
       write(*,'(1x,''|-----------------------------------------------|'')')
       write(*,'(1x,''|       [s]ave   <[c]ontinue>   [r]eplot        |'')') 
       write(*,'(1x,''|_______________________________________________|'')')
       write(*,'(/1x,'' enter option number > '',$)')
       read(*,'(a1)')ans
       if((ans.eq.'S').or.(ans.eq.'s'))then
         dev='psc';ok=.false.;filename='fig11.ps'
	  call plsdev(dev)
         call plsfnam(filename)
         goto 10
       elseif((ans.eq.'R').or.(ans.eq.'r'))then
         ok=.true.
         goto 10
       elseif((ans.eq.'C').or.(ans.eq.'c').or.(ans.eq.' ').or.
     *                             (ans.eq.'Y').or.(ans.eq.'y'))then
         return
	endif
c
	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine plt3(mxy,mxs,indx,xrc,f,idate,flag,nmax,sernum,n,ident,string,plotall)
c
c	argument
	integer mxy,mxs,nmax,sernum,n
	integer idate(mxs,3)
	real*4 indx(mxy),flag,f(mxy),xrc(mxy)
	character ident(mxs)*8
	character string*30
	logical plotall
c	local
	integer p,j,k, minx, maxx
	integer red,blue,green
	real*4 miny, maxy
	real*4 wrk1(5000),wrk2(5000),wrk3(5000),wrk4(5000)
	character dev*80
	logical ok,bg
c
       dev='Mac1';ok=.true.
       red=255;green=255;blue=255
       bg=.false.
c
10	if(ok)call plsdev(dev)
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	
	miny=9999.;maxy=-9999.;minx=9999;maxx=-9999
	do p=1,nmax
	   if(xrc(p).ne.flag.and.xrc(p).gt.maxy)maxy=xrc(p)
	   if(f(p).ne.flag.and.f(p).gt.maxy)maxy=f(p)
	   wrk1(p)=real(p)
	enddo
c
c initialize plplot
	call plinit()
c set lable color
	call plcol(15)   !color box and lables
c define the viewport one column by one rows
	call plssub(1,2)
	call plenv(real(1),real(nmax),0.,maxy,0,2)
	call pllab('Biological Year','Y',string)
	call plcol(1)   !color box and lables
	call plline(n,wrk1,xrc)
	call plcol(9)   !color box and lables
	call plwid(2)
	call plline(nmax,wrk1,f)
	call plwid(1)
c
c bottom plot
	miny=9999.; maxy=-9999.; minx=9999; maxx=-9999; j=idate(sernum,1)
	do p=1,n
	      if(indx(p).ne.flag.and.indx(p).gt.maxy)maxy=indx(p)
	      if(indx(p).ne.flag.and.indx(p).lt.miny)miny=indx(p)
	      wrk1(p)=real(j)
	      j=j+1
	enddo
c
c set lable color
	call plcol(15)   !color box and lables
	call plenv(real(idate(sernum,1)),real(idate(sernum,2)),miny,maxy,0,2)
	call pllab('Calendar Year','Y',ident(sernum))
	call plcol(1)   !color box and lables
	call plline(n,wrk1,indx)
	call plend()
c	
       ans='C '
	dev='Mac1'
       write(*,'(/1x,''|************* plot options menu ***************|'')')
       write(*,'(1x,''|-----------------------------------------------|'')')
       write(*,'(1x,''| [s]ave  [c]ontinue  [r]eplot  [n]o more plots |'')') 
       write(*,'(1x,''|_______________________________________________|'')')
       write(*,'(/1x,'' enter option number > '',$)')
       read(*,'(a1)')ans
       if((ans.eq.'S').or.(ans.eq.'s'))then
         dev='psc';ok=.false.;filename='fig11.ps'
	  call plsdev(dev)
         call plsfnam(filename)
         goto 10
       elseif((ans.eq.'R').or.(ans.eq.'r'))then
         ok=.true.
         goto 10
       elseif((ans.eq.'C').or.(ans.eq.'c').or.(ans.eq.' ').or.
     *                             (ans.eq.'Y').or.(ans.eq.'y'))then
         return
       elseif((ans.eq.'N').or.(ans.eq.'n'))then
         plotall=.false.
         return
	endif
c
	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine plt4(mxlen,maxit,rcscrv,nmax,rcnt)
c plot all regional curves and their solutions
	integer mxlen,maxit,rcnt
	real*4 rcscrv(mxlen,2,maxit)
c	local
	integer i,p,j,k,nrows,ncols,nplts
	integer red,blue,green
	real*4 miny, maxy, minx, maxx
	real*4 x(5000),y1(5000),y2(5000)
	character string*16,dev*3
	logical ok,bg
c
       dev='Mac1';ok=.true.
       red=255;green=255;blue=255
       bg=.false.
c
10	if(ok)call plsdev(dev)
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	
	miny=9999.;maxy=-9999.
	minx=0.;maxx=float(nmax)
	do k=1,rcnt
	  do j=1,2
	    do p=1,mxlen
	      if(rcscrv(p,j,k).gt.maxy)maxy=rcscrv(p,j,k)
	      if(rcscrv(p,j,k).lt.miny.and.rcscrv(p,j,k).ne.0.0)miny=rcscrv(p,j,k)
	    enddo
	  enddo
	enddo
c
	if(rcnt.le.8)then
	   nrows=2;ncols=4
	elseif(rcnt.gt.8.and.rcnt.le.10)then
	   nrows=2;ncols=5
	elseif(rcnt.gt.10.and.rcnt.le.15)then
	   nrows=3;ncols=5
	elseif(rnct.gt.15.and.rcnt.le.20)then
	   nrows=4;ncols=5
	else
	   write(*,'(/1x,''Too many Regional curves to plot. Try '')')
	   write(*,'(1x,''plotting them outside this program.'')')
	   write(*,'(1x,''Hit RETURN to Continue'')')
	   pause
	   return
	endif
c	
	string='Iteration #     '
	nplts=nrows*ncols
c
c initialize plplot
	call plinit()
	call plssub(nrows,ncols)
	do i=1,rcnt
	    	  do p=1,nmax
	    	       y1(p)=rcscrv(p,1,i)
	              y2(p)=rcscrv(p,2,i)
	              if(i.eq.1)x(p)=real(p)
	         enddo
	         call plcol(15)   !color box and lables
	         call plenv(minx,maxx,miny,maxy,0,2)
	         if(i.eq.1)then
	             call pllab('Biological Year','Y','Initial RC and Fit')
	         else
	            write(string(14:16),'(i3)')i-1
	            call pllab('Biological Year','Y',string)
	         endif
	         call plcol(9)
	         call plwid(1)	     
	         call plline(nmax,x,y1)
	         call plcol(1)   !color box and lables
	         call plwid(2)	     
	         call plline(nmax,x,y2)
	         call plwid(1)	
	enddo
	if(rcnt.lt.nplts)then
	    do i=1,(nplts-rcnt)
	       call pladv()
	    enddo
	endif
c
20	call plend()	    	    
c	
       ans='C '
	dev='Mac1'
       write(*,'(/1x,''|************* plot options menu ***************|'')')
       write(*,'(1x,''|-----------------------------------------------|'')')
       write(*,'(1x,''|       [s]ave   <[c]ontinue>   [r]eplot        |'')') 
       write(*,'(1x,''|_______________________________________________|'')')
       write(*,'(/1x,'' enter option number > '',$)')
       read(*,'(a1)')ans
       if((ans.eq.'S').or.(ans.eq.'s'))then
         dev='psc';ok=.false.;filename='fig11.ps'
	  call plsdev(dev)
         call plsfnam(filename)
         goto 10
       elseif((ans.eq.'R').or.(ans.eq.'r'))then
         ok=.true.
         goto 10
       elseif((ans.eq.'C').or.(ans.eq.'c').or.(ans.eq.' ').or.
     *                             (ans.eq.'Y').or.(ans.eq.'y'))then
         return
	endif
c
	return
	end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine slpt(mxy,mxs,x,idate,nc)
c performs adaptive power transformation on matrix x alla cook and peters 1997, holocene
c pjk 2011, SU
c argument decs
	integer mxy,mxs,nc
	real*4 x(mxy,mxs)
	integer idate(mxs,3)
	
c local declarations
	integer p,j,k
	integer maxdat,maxser
	parameter(maxdat=5000,maxser=3000)!dimension must equal mxy
	integer num(maxser)
	real*4 mr(maxdat,maxser),sr(maxdat,maxser) 
	real*4 sl(maxser),yi(maxser)
c	
	do k=1,nc
	       num(k)=(idate(k,2)-idate(k,1))+1
		do p=2,num(k)
	  		mr(p-1,k)=log( ((x(p,k)+x(p-1,k))/2.) )
	  		if(abs(x(p,k)-x(p-1,k)).eq.0)then
	  		    sr(p-1,k)=log(0.001)
	  		else
	  		    sr(p-1,k)=log( abs(x(p,k)-x(p-1,k))   ) 
	  		endif
c	  		write(*,*)k,p,mr(p-1,k),sr(p-1,k)
		enddo
		call trend2(num(k)-1,mr(1,k),sr(1,k),sl(k),yi(k))
		pt=1.0-sl(k)
		do p=1,num(k)
	   		x(p,k)=x(p,k)**pt
		enddo
	enddo
	call scatplt(maxdat,maxser,mr,sr,sl,yi,num,nc)
	return
	end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine scatplt(mxy,mxs,x1,y1,slope,int,num,nc)
c create a scatter plot of x vs y and thier ols solution
	integer mxy,mxs,num(mxs),nc
	real*4 x1(mxy,mxs),y1(mxy,mxs),slope(mxs),int(mxs)
c	local
	integer p,j,k
	integer red,blue,green
	real*4 miny, maxy, minx, maxx
	real*4 x2(5000),y2(5000),y3(5000)
	character string*16,dev*3,filename
	logical ok,bg
c
       dev='Mac1';ok=.true.
       red=255;green=255;blue=255
       bg=.false.
c
10	if(ok)call plsdev(dev)
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	
	miny=9999.;maxy=-9999.
	minx=miny ;maxx=maxy
	do k=1,nc
	   do p=1,num(k)-1
	      if(y1(p,k).gt.maxy)maxy=y1(p,k)
	      if(y1(p,k).lt.miny)miny=y1(p,k)
	      if(x1(p,k).gt.maxx)maxx=x1(p,k)
	      if(x1(p,k).lt.minx)minx=x1(p,k)
	  enddo
	enddo
c
c initialize plplot
	call plinit()
	call plcol(15)   !color box and lables
	call plenv(minx,maxx,miny,maxy,0,2)
	do k=1,nc
	    call plcol(9)
	    do p=1,num(k)-1
	       x2(p)=x1(p,k)
	       y2(p)=y1(p,k)
	    enddo
	    call plpoin(num(k)-1,x2,y2,1)
	    call plcol(1)
	    do p=1,num(k)-1
	      y3(p)=x2(p)*slope(k)+int(k)
	    enddo
	    call plline(num(k)-1,x2,y3)
	enddo
20	call plend()	    	    
c	
       ans='C '
	dev='Mac1'
       write(*,'(/1x,''|************* plot options menu ***************|'')')
       write(*,'(1x,''|-----------------------------------------------|'')')
       write(*,'(1x,''|       [s]ave   <[c]ontinue>   [r]eplot        |'')') 
       write(*,'(1x,''|_______________________________________________|'')')
       write(*,'(/1x,'' enter option number > '',$)')
       read(*,'(a1)')ans
       if((ans.eq.'S').or.(ans.eq.'s'))then
         dev='psc';ok=.false.;filename='fig1b.ps'
	  call plsdev(dev)
         call plsfnam(filename)
         goto 10
       elseif((ans.eq.'R').or.(ans.eq.'r'))then
         ok=.true.
         goto 10
       elseif((ans.eq.'C').or.(ans.eq.'c').or.(ans.eq.' ').or.
     *                             (ans.eq.'Y').or.(ans.eq.'y'))then
         return
	endif
c
	return
	end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*rtn trend2: fits linear regression line to data pairs
*
* programmed (c) may 1987 from trend by
* richard l holmes, holzbiologie - bfh, universitaet hamburg - germany
* modified (c) 04 aug 1990
*
*  n:  number of values in each array
*  x:  array of horizontal coordinates
*  y:  array of vertical coordinates
*  sl: slope of regression
*  yi: y-axis intercept
*
*  equation:  y = sl * x + yi
*
      subroutine trend2(n,x,y,sl,yi)
      real x(1),y(1)
      an=float(n)
      t=0.
      w=0.
      v=0.
      z=0.
      do 10 i=1,n
      z=z+x(i)
      t=t+x(i)*x(i)
      v=v+y(i)*x(i)
   10 w=w+y(i)
      ymean=w/an
      if((t-z*z/an) .ne. 0.)then
      sl=(v-w*z/an)/(t-z*z/an)
      else
      sl=0.
      endif
      yi=ymean-sl*z/an
      end ! trend2
