	program SignalFree
c
	integer iu,ou,mxd,mxr
	parameter(mxd=10000,mxr=10000,iu=69,ou=68)
	integer p,j,k,ifl,ls,nrec
	integer mniyr,mxlyr,mnage,mxage,ici,lci
	integer age(mxr),iyr(mxr),lyr(mxr),cyr1,cyr2,numsft
	real   rrw(mxr,mxd),indx(mxr,mxd),crn(mxd),sft(mxr,mxd)
	real   y1(mxd),y2(mxd),nc(mxd),nrrw(mxd),sftcrn(mxd)
	real   sum(mxd),sumsft,sumrrw
	character*1 z
	character*8 id(mxr)
	character*40 infl,outfl
c
	write(*,'(1x,''Enter input file name : '',$)')
	read(*,'(a)')infl
	do p=1,40-1
	  read(infl(p:p),'(a)')z
	  if(z.eq.' ')goto 5
	end do
5	outfl=infl(1:p-1)//'_out'
	write(*,'(/1x,''input file  = '',a)')infl
	write(*,'(1x,''output file = '',a)')outfl
c
	mniyr=mxd
	mxlyr=-1.0*mniyr
	ls=100
	do p=1,mxd
	  sum(p)=0.00
	  nrrw(p)=0.00
	end do
c open and read
	open(unit=iu,file=infl,status='old',readonly)
	p=1
10	call trrw(mxd,id(p),iyr(p),age(p),rrw(p,1),nc,ifl,iu,iuo,'MR') 
c	  write(*,'(1x,''ifl= '',i5)')ifl
	  if(ifl.gt.0)then 
	    nrec=ifl
           lyr(p)=iyr(p)+age(p)-1
           k=1
	    do j=iyr(p),lyr(p)
	      nrrw(j)=nrrw(j)+1.0
	      k=k+1
	    end do	    
	    if(iyr(p).lt.mniyr)mniyr=iyr(p)
	    if(lyr(p).gt.mxlyr)mxlyr=lyr(p)
	    if(age(p).lt.mnage)mnage=age(p)
	    if(age(p).gt.mxage)mxage=age(p)
	    p=p+1
	    goto 10
	  else
	    write(*,'(/1x,''INPUT file read without errors'')')
	    p=p-1
	  endif
c
	write(*,'(/1x,'' Number of series read = '',i5)')nrec
	write(*,'(1x,''series   id        iyr   lyr   age   '')')
	do p=1,nrec
	   write(*,'(1x,i5,3x,a8,3i6)')p,id(p),iyr(p),lyr(p),age(p)
 	end do
 	close(iu)
c compute common period
	call comint(mxr,iyr,id,age,nrec,ici,lci)
c spline all rrw series and make chronolonogy
	do p=1,nrec
	   k=1
	   do j=iyr(p),lyr(p)
	      y1(k)=rrw(p,j)
	      nrrw(j)=nrrw(j)+1
	      k=k+1
	   end do
	   call spline(mxd,age(p),y1,y2,ls)
	   k=1
	   do j=iyr(p),lyr(p)
	      indx(p,j)=y1(k)/y2(k)
	      sum(j)=sum(j)+indx(p,j)  !sum indices
	      k=k+1
	   end do
	end do
	do p=mniyr,mxlyr
	   crn(p)=sum(p)/nrrw(p)
	end do
	
c determine where nrrw.ge.2 
	do p=mniyr,mxlyr
	  if(nrrw(p).gt.1)then
	    cyr1=p
	    goto 20
	  end if
	end do
20	do p=mxlyr,mniyr,-1
	  if(nrrw(p).gt.1)then
	    cyr2=p
	    goto 30
	  end if
	end do
30	continue
	numsft=cyr2-cyr1+1
	do k=1,mxd
	  sum(k)=0.00
	  nrrw(k)=0.00
	end do
c
	do p=1,nrec
	    do k=cyr1,cyr2
	       sft(p,k)=rrw(p,k)/crn(k)
	    end do
	    sumrrw=asum(rrw(p,1),cyr1,cyr2)
	    sumsft=asum(sft(p,1),cyr1,cyr2)
	    do k=cyr1,cyr2
	       sft(p,k)=sft(p,k)*sumrrw/sumsft
	       nrrw(k)=nrrw(k)+1.0
	    end do
c spline the new sft array	    
	    j=1
	    do k=cyr1,cyr2
	       y1(j)=sft(p,k)
	       j=j+1
	    end do
	    call spline(mxd,numsft,y1,y2,ls)
	    j=1
	    do k=cyr1,cyr2
	         sft(p,k)=y1(j)/y2(j)
	         sum(k)=sum(k)+sft(p,k)
	         j=j+1
	    end do
	end do
c make new sft chronology and signal free chronology
	do k=cyr1,cyr2
	    sftcrn(k)=sum(k)/nrrw(k)
	    crn(k)=crn(k)*sftcrn(k)
	end do
c	  
	stop
	end
