! Copyright (C) 2013, Thomas M. Melvin and Keith R. Briffa, see 
! the GNU General Public License.
      MODULE  rcs    
      USE crustutil
      IMPLICIT NONE
      CHARACTER(30),DIMENSION(fin) :: rcnam
      CONTAINS   
!--------------------------------------------------------
      SUBROUTINE rcs_val()  
      IMPLICIT NONE                 
      rcnam(1:30)= &
       (/"Exit RCS menu                 ", &
         "Save this plot                ", &
         "( 3)                          ", &
         "( 4) Fig1 Torn All v Living   ", &
         "( 5) Fig2 Long Step RCS       ", &
         "( 6) Fig3 Long Step CRN       ", &
         "( 7) Fig4 Short Step RCS      ", &
         "( 8) Fig5 Short Step CRN      ", &
         "( 9) Fig6 HB Rand RCS         ", &
         "(10) Fig7 HB Rand CRN         ", &
         "(11) Fig8 Short Slope CRN     ", &
         "(12) Fig9 Ages and Steps      ", &
         "(13)                          ", &
         "(14)                          ", &
         "(15)                          ", &
         "(16)                          ", &
         "(17)                          ", &
         "(18)                          ", &
         "(19)                          ", &
         "(20)                          ", &
         "(21)                          ", &
         "(22)                          ", &
         "(23)                          ", &
         "(24)                          ", &
         "(25)                          ", &
         "(26)                          ", &
         "(27)                          ", &
         "(28)                          ", &
         "(29)                          ", &
         "(30)                          "/)
      rcnam(31:60)= &
       (/"(31) SM01 Multi-Slope         ", &
         "(32) SM02 Short Slope RCS     ", &
         "(33) SM03 RCS slop no mean    ", &
         "(34) SM04 CRN slop no mean    ", &
         "(35) SM05 RCS mean no slop    ", &
         "(36) SM06 CRN mean no slop    ", &
         "(37) SM07 150 year            ", &
         "(38) SM08 200 year            ", &
         "(39) SM09 250 year            ", &
         "(40) SM10 300 year            ", &
         "(41) SM11 350 year            ", &
         "(42) SM12 Ages and Slopes     ", &
         "(43)                          ", &
         "(44)                          ", &
         "(45)                          ", &
         "(46)                          ", &
         "(47)                          ", &
         "(48)                          ", &
         "(49)                          ", &
         "(50)                          ", &
         "(51)                          ", &
         "(52)                          ", &
         "(53)                          ", &
         "(54)                          ", &
         "(55)                          ", &
         "(56)                          ", &
         "(57)                          ", &
         "(58)                          ", &
         "(59)                          ", & 
         "(60)                          "/) 
      rcnam(61:90)= &
       (/"(61) 2RC1 Pith Offset         ", &
         "(62) 2RC2 SDev Indices        ", &
         "(63) 2RC3 U Curves            ", &
         "(64) 2RC4 RCS EPS             ", &
         "(65) 2RC5 1, 2 and 3 RCS      ", &
         "(66) 2RC6 Yam Sl/Mn CRN       ", &
         "(67) 2RC7 China Sl/Mn CRN     ", &
         "(68) 2RC8 MRCS Random         ", &
         "(69)                          ", &
         "(70)                          ", &
         "(71) R2SM1 Yanmal SD PTrans   ", &
         "(72) R2SM2 Yamal Sl/Mn SDev   ", &
         "(73) R2SM3 China Sl/Mn SDev   ", &
         "(74)                          ", &
         "(75)                          ", & 
         "(76)                          ", &
         "(77)                          ", &
         "(78)                          ", &
         "(79)                          ", &
         "(80)                          ", &
         "(81)                          ", &
         "(82)                          ", &
         "(83)                          ", &
         "(84)                          ", &
         "(85)                          ", &
         "(86)                          ", &
         "(87)                          ", &
         "(88)                          ", &
         "(89)                          ", &
         "(90)                          "/)
      RETURN 
      END SUBROUTINE rcs_val
!--------------------------------------------------------
      SUBROUTINE rcsv(ref1,plot)  
      IMPLICIT NONE                 
      INTEGER,INTENT(IN) :: ref1   ! Chosen action
      INTEGER,INTENT(IN) :: plot   ! Plot number
      INTEGER            :: i
      CALL ERASE()
      SELECT CASE (ref1)
      CASE (2)    ! Save plot as .ps file
        figm="rcs/rcs0"
        portrait=plot.GE.61.AND.plot.LE.73  ! Portrait 
        CALL open_ps(plot,9)
        SELECT CASE (plot)
          CASE ( 4) ; CALL R1Fig01d() 
          CASE ( 5) ; CALL R1Fig02d(3)
          CASE ( 6) ; CALL R1Fig03d(3)
          CASE ( 7) ; CALL R1Fig02d(1)
          CASE ( 8) ; CALL R1Fig03d(1)
          CASE ( 9) ; CALL R1Fig08d()
          CASE (10) ; CALL R1Fig07d()
          CASE (11) ; CALL R1Fig03d(2)
          CASE (12) ; CALL R1Fig09d(1)     
          CASE (31) ; CALL R1SM01d()
          CASE (32) ; CALL R1Fig02d(2)
          CASE (33) ; CALL R1Fig02d(4)
          CASE (34) ; CALL R1Fig03d(4)
          CASE (35) ; CALL R1Fig02d(5)
          CASE (36) ; CALL R1Fig03d(5)
          CASE (37:41) ; CALL R1Fig03d(6)
          CASE (42) ; CALL R1Fig09d(2)     
          CASE (61) ; CALL R2Fig1d() 
          CASE (62) ; CALL R2Fig2d()
          CASE (63) ; CALL R2Fig3d()
          CASE (64) ; CALL R2Fig4d()
          CASE (65) ; CALL R2Fig5d()
          CASE (66) ; CALL smerrd(1)     
          CASE (67) ; CALL smerrd(3)     
          CASE (68) ; CALL R2Fig8d()
          CASE (71) ; CALL R2SM1d() 
          CASE (72) ; CALL smerrd(2)     
          CASE (73) ; CALL smerrd(4)     
        END SELECT
        CALL plot_psend()
      CASE ( 4) ; CALL R1Fig01()  ; CALL R1Fig01d()   ! Yamal All v Living 
      CASE ( 5) ; CALL RCFig02()  ; CALL R1Fig02d(3)  ! Step on long RCS
      CASE ( 6) ; CALL RCFig02()  ; CALL R1Fig03d(3)  ! Step on long CRN
      CASE ( 7) ; CALL R1Fig03()  ; CALL R1Fig02d(1)  ! Step on short RCS
      CASE ( 8) ; CALL R1Fig03()  ; CALL R1Fig03d(1)  ! Step on short CRN
      CASE ( 9) ; CALL R1Fig07()  ; CALL R1Fig08d()   ! Slope long/short RCS
      CASE (10) ; CALL R1Fig07()  ; CALL R1Fig07d()   ! Slope long/short CRN
      CASE (11) ; CALL RCFig05(1) ; CALL R1Fig03d(2)  ! Slope on short CRN
      CASE (12) ; CALL R1Fig09(1) ; CALL R1Fig09d(1)  ! Step on short CRN
      CASE (31) ; CALL R1SM01()   ; CALL R1SM01d()    ! Sig-free and slopes 
      CASE (32) ; CALL RCFig05(1) ; CALL R1Fig02d(2)  ! Slope on short RCS
      CASE (33) ; CALL RCFig05(2) ; CALL R1Fig02d(4)  ! Slope no mean RCS
      CASE (34) ; CALL RCFig05(2) ; CALL R1Fig03d(4)  ! Slope no mean CRN
      CASE (35) ; CALL RCFig05(3) ; CALL R1Fig02d(5)  ! Mean no slope RCS
      CASE (36) ; CALL RCFig05(3) ; CALL R1Fig03d(5)  ! Mean no slope CRN
      CASE (37) ; CALL RC_sm21(2) ; CALL R1Fig03d(6)  ! Step on short CRN
      CASE (38) ; CALL RC_sm21(3) ; CALL R1Fig03d(6)  ! Step on short CRN
      CASE (39) ; CALL RC_sm21(4) ; CALL R1Fig03d(6)  ! Step on short CRN
      CASE (40) ; CALL RC_sm21(5) ; CALL R1Fig03d(6)  ! Step on short CRN
      CASE (41) ; CALL RC_sm21(6) ; CALL R1Fig03d(6)  ! Step on short CRN
      CASE (42) ; CALL R1Fig09(2) ; CALL R1Fig09d(2)  ! Slope on short CRN
      CASE (61) ; CALL R2Fig1() ; CALL R2Fig1d()    ! RCS Pith Offset 
      CASE (62) ; CALL R2Fig2() ; CALL R2Fig2d()    ! RC-ratio CRN and SDev
      CASE (63) ; CALL R2Fig3() ; CALL R2Fig3d()    ! RC-index distribution
      CASE (64) ; CALL R2Fig4() ; CALL R2Fig4d()    ! EPSS figure
      CASE (65) ; CALL R2Fig5() ; CALL R2Fig5d()    ! Random RCS
      CASE (66) ; CALL smerr(1) ; CALL smerrd(1)    ! Slope and mean error
      CASE (67) ; CALL smerr(3) ; CALL smerrd(3)    ! Slope and mean error
      CASE (68) ; CALL R2Fig8() ; CALL R2Fig8d()    ! Random RCS
      CASE (71) ; CALL R2SM1()  ; CALL R2SM1d()     ! Yamal PT SDev
      CASE (72) ; CALL smerr(2) ; CALL smerrd(2)    ! Slope and mean error
      CASE (73) ; CALL smerr(4) ; CALL smerrd(4)    ! Slope and mean error
      ENDSELECT
      b(272)%ok=TR
      IDO: DO i=1,50
        b(271)%on=TR ; CALL but_draw(272,"")
        CALL mouse_click(3,272,272) ; IF (mous.EQ.272) EXIT IDO
      ENDDO IDO
      b(272)%ok=FA
      RETURN 
      END SUBROUTINE rcsv
!----------------------------------------------------------------
      SUBROUTINE R1Fig03()  ! Reduced rings final years 
      IMPLICIT NONE
      INTEGER              :: i,j,p,q,r
      REAL(8),DIMENSION(3) :: fac
      OPEN(74,FILE="zzz.prn",IOSTAT=ios,STATUS="REPLACE")
      CALL det_default() ; idt=-2
      fac=(/1.D0,0.7D0,0.4D0/) ; cf=1 
      DO j=1,3
        nc=0 ; CALL read_rft("../../raw/rcs/finnmrg.raw")
        DO i=1,nc
          p=MAX(fy(i),1920) ; q=ly(i)
          IF (p.LT.q) THEN
            p=ad(i)+(p-fy(i)) ; q=ad(i+1)-1
            x(p:q)=x(p:q)*fac(j)
          ENDIF
        ENDDO
        IF (j.EQ.2) THEN
          cf=cf+5 ; idt=0 ; CALL detrend() ; cf=cf-5  ! cf=8
        ELSEIF (j.EQ.3) THEN
          cf=cf+4 ; idt=0 ; CALL detrend() ; cf=cf-4  ! cf=9
        ELSE
          idt=0 ; CALL detrend()         ! cf=1
        ENDIF
        idt=-2 ; cf=cf+1 ; sfo=1 ; CALL detrend()
        r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
        num(1:r,cf+10)=mcnt(1:r,mx)
        cf=cf+1 ; sfo=2 ; CALL detrend()
        r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
        num(1:r,cf+10)=mcnt(1:r,mx)
        WRITE(74,'(4I6)') j,cfy(cf-2),cly(cf-2),cyr(cf-2)
        WRITE(74,'(4I6)') cf,cfy(cf-1),cly(cf-1),cyr(cf-1)
        WRITE(74,'(4I6)') cf,cfy(cf),cly(cf),cyr(cf)
      ENDDO
      r=cyr(1)
      crn(1:r,36)=crn(1:r,6)/crn(1:r,4)
      crn(1:r,37)=crn(1:r,2)/crn(1:r,4)
      crn(1:r,38)=crn(1:r,7)/crn(1:r,5)
      crn(1:r,39)=crn(1:r,3)/crn(1:r,5)
      DO i=1,12
        WRITE(74,'(4I6)') i,cfy(i),cly(i),cyr(i)
      ENDDO
      CLOSE(74)
      cnam(30)="1920+ Rings *0.4"
      cnam(31)="Measured Values"
      cnam(32)="1920+ Rings *0.7"
      cnam(33)="a) Finnish Trees - RCS Curves"
      cnam(34)="a) Finnish Trees - Ring-Width"
      RETURN
      END SUBROUTINE R1Fig03
!--------------------------------------------------------------
      SUBROUTINE RCFig05(ref1)  ! Slope adjust no change to mean
      IMPLICIT NONE
      INTEGER,INTENT(IN)   :: ref1 
      INTEGER              :: i,j,p,q,r,u,v
      REAL(8)              :: ra,rb
      REAL(8),DIMENSION(3) :: fac
      CALL det_default() ; idt=-2
      OPEN(74,FILE="zzz.prn",IOSTAT=ios,STATUS="REPLACE")
      fac=(/+0.001D0,0.D0,-0.001D0/) ; cf=1 
      DO j=1,3
        nc=0 ; CALL read_rft("../../raw/rcs/finnmrg.raw")
        CALL det_crnfy() ; r=cyr(1)
        crn(1:r,49+j)=(/(1.D0+DBLE(i)*fac(j),i=-r/2,r-r/2-1)/)
        DO i=1,nc
          r=yr(i) ; p=ad(i) ; q=p+r-1 
          u=fy(i)-cfy(1)+1 ; v=u+r-1 ; ra=SUM(x(p:q))
          IF (ref1.EQ.2) THEN       ! Slope no mean
            x(p:q)=x(p:q)*crn(u:v,49+j) ; rb=SUM(x(p:q))
            x(p:q)=x(p:q)*ra/rb
          ELSEIF (ref1.EQ.3) THEN   ! Mean no slope
            rb=SUM(x(p:q)*crn(u:v,49+j))
            x(p:q)=x(p:q)*ra/rb
          ELSEIF (ref1.EQ.1) THEN   ! Slope and mean
            x(p:q)=x(p:q)*crn(u:v,49+j) 
          ENDIF
        ENDDO
        IF (j.EQ.2) THEN
          cf=cf+5 ; idt=0 ; CALL detrend() ; cf=cf-5  ! cf=8
        ELSEIF (j.EQ.3) THEN
          cf=cf+4 ; idt=0 ; CALL detrend() ; cf=cf-4  ! cf=9
        ELSE
          idt=0 ; CALL detrend()         ! cf=1
        ENDIF
        idt=-2 ; cf=cf+1 ; sfo=1 ; CALL detrend()
        r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
        num(1:r,cf+10)=mcnt(1:r,mx)
        cf=cf+1 ; sfo=2 ; CALL detrend()
        r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
        num(1:r,cf+10)=mcnt(1:r,mx)
        WRITE(74,'(4I6)') j,cfy(cf-2),cly(cf-2),cyr(cf-2)
        WRITE(74,'(4I6)') cf,cfy(cf-1),cly(cf-1),cyr(cf-1)
        WRITE(74,'(4I6)') cf,cfy(cf),cly(cf),cyr(cf)
      ENDDO
      r=cyr(1)  ! Count weighted centre is year 1859.5 
      crn(1:r,50)=crn(1:r,50)/(SUM(crn(310:311,50))/2.D0)   
      crn(1:r,52)=crn(1:r,52)/(SUM(crn(310:311,52))/2.D0)  
      DO i=1,12
        WRITE(74,'(4I6)') i,cfy(i),cly(i),cyr(i)
      ENDDO
      DO i=1,r
        WRITE(74,'(3I7)') cfy(1)-1+i,i,SUM(num(1:i,1))
      ENDDO  
      CLOSE(74)
      crn(1:r,36)=crn(1:r,6)/crn(1:r,4)
      crn(1:r,37)=crn(1:r,2)/crn(1:r,4)
      crn(1:r,38)=crn(1:r,7)/crn(1:r,5)
      crn(1:r,39)=crn(1:r,3)/crn(1:r,5)
      CLOSE(74)
      cnam(30)="Slope -0.001"
      cnam(31)="Slope +0.001"
      cnam(32)="Measured Values"
      IF (ref1.EQ.2) THEN
        cnam(33)="a) No mean Finnish Trees - RCS Curves"
        cnam(34)="a) No mean Finnish Trees - Ring-Width"
      ELSEIF (ref1.EQ.3) THEN
        cnam(33)="a) No slope Finnish Trees - RCS Curves"
        cnam(34)="a) No slope Finnish Trees - Ring-Width"
      ELSEIF (ref1.EQ.1) THEN
        cnam(33)="a) Finnish Trees - RCS Curves"
        cnam(34)="a) Finnish Trees - Ring-Width"
      ENDIF
      RETURN
      END SUBROUTINE RCFig05
!--------------------------------------------------------------
      SUBROUTINE RCFig02()  ! Reduced rings final years 
      IMPLICIT NONE
      INTEGER              :: i,j,p,q,r
      REAL(8),DIMENSION(3) :: fac
      OPEN(74,FILE="zzz.prn",IOSTAT=ios,STATUS="REPLACE")
      CALL det_default() ; idt=-2
      fac=(/1.4D0,1.D0,0.6D0/) ; cf=1 
      DO j=1,3
        nc=0 ; CALL read_rft("../../raw/rcs/yml-all.raw")
        DO i=1,nc              ! For each tree
          p=MAX(fy(i),1880) ; q=ly(i)
          IF (p.LT.q) THEN     ! Adjust after 1880
            p=ad(i)+(p-fy(i)) ; q=ad(i+1)-1
            x(p:q)=x(p:q)*fac(j)
          ENDIF
        ENDDO
        IF (j.EQ.2) THEN
          cf=cf+5 ; idt=0 ; CALL detrend() ; cf=cf-5  ! cf=8
        ELSEIF (j.EQ.3) THEN
          cf=cf+4 ; idt=0 ; CALL detrend() ; cf=cf-4  ! cf=9
        ELSE
          idt=0 ; CALL detrend()         ! cf=1
        ENDIF
        idt=-2 ; cf=cf+1 ; sfo=1 ; CALL detrend()
        r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
        num(1:r,cf+10)=mcnt(1:r,mx)
        cf=cf+1 ; sfo=2 ; CALL detrend()
        r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
        num(1:r,cf+10)=mcnt(1:r,mx)
        WRITE(74,'(4I6)') j,cfy(cf-2),cly(cf-2),cyr(cf-2)
        WRITE(74,'(4I6)') cf,cfy(cf-1),cly(cf-1),cyr(cf-1)
        WRITE(74,'(4I6)') cf,cfy(cf),cly(cf),cyr(cf)
      ENDDO
      DO i=1,12
        WRITE(74,'(4I6)') i,cfy(i),cly(i),cyr(i)
      ENDDO
      r=cyr(1)
      crn(1:r,36)=crn(1:r,6)/crn(1:r,4)
      crn(1:r,37)=crn(1:r,2)/crn(1:r,4)
      crn(1:r,38)=crn(1:r,7)/crn(1:r,5)
      crn(1:r,39)=crn(1:r,3)/crn(1:r,5)
      DO i=1,9
        CALL splinet(r,crn(1:r,i),50,crn(1:r,i)) 
      ENDDO
      CLOSE(74)
      cnam(30)="1880+ Rings *0.6"
      cnam(31)="1880+ Rings *1.4"
      cnam(32)="Measured Values"
      cnam(33)="a) Yamal Trees - RCS Curves"
      cnam(34)="a) Yamal Trees - Ring-Width"        
      RETURN
      END SUBROUTINE RCFig02
!--------------------------------------------------------------
      SUBROUTINE R1Fig03d(ref1)  
      IMPLICIT NONE                 
      INTEGER,INTENT(IN)     :: ref1
      REAL(8),DIMENSION(mxy) :: wka
      REAL(8)                :: ra,rb,rc,rd,re,rf  
      CHARACTER(12)          :: jnam
      INTEGER                :: i=1,r,p,q,y1,y2
      IF     (ref1.EQ.1) THEN
        ra=0.D0   ; rb=1.2D0  ; rc=0.1D0 ; rd=1.7D0
        re=0.6D0  ; rf=1.3D0
        jnam="R1fig5.col" ; wnam(49)="Finnish trees"
      ELSEIF (ref1.EQ.2.OR.ref1.EQ.4.OR.ref1.EQ.5) THEN
        ra=0.1D0  ; rb=1.2D0  ; rc=0.2D0 ; rd=1.99D0
        re=0.65D0  ; rf=1.35D0 ; wnam(49)="Finnish trees"
        IF (ref1.EQ.2) jnam="R1fig8.col"
        IF (ref1.EQ.4) jnam="R1sm4.col"
        IF (ref1.EQ.5) jnam="R1sm6.col"
      ELSEIF (ref1.EQ.3) THEN
        ra=0.25D0 ; rb=1.5D0  ; rc=0.3D0 ; rd=2.4D0
        re=0.64D0  ; rf=1.32D0 
        jnam="R1fig3.col" ; wnam(49)="Yamal trees"
      ELSEIF (ref1.EQ.6) THEN
        ra=0.1D0 ; rb=1.6D0  ; rc=0.2D0 ; rd=1.8D0
        re=0.6D0  ; rf=1.3D0 
      ENDIF
      y1=cfy(1) ; y2=cly(1)
      wka(1:cyr(1))=(/(DBLE(i),i=cfy(1),cly(1))/)   
      p=y1-cfy(1)+1 ; q=y2-cfy(1)+1 ; r=q-p+1 
      grl=200 ; grr=2000 ; grt=140 ; grb=430
      CALL TICKS(5,'X')            ! No X ticks
      CALL plot_trees(r,num(p:q,1))  
      CALL LABELS('NONE','X')
      CALL NAME('','X')            ! Axis name
      CALL NAME('Index Value','Y') ! Axis name
      CALL NAME('Ring Width','Y')  ! Axis name
      CALL tombox(y1,y2,ra,rb)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(blue)  ; CALL CURVE(wka(p:q),crn(p:q,9),r)
      CALL MESSAG(cnam(30),grl+1400,grt-40)
      CALL SETCLR(red)   ; CALL CURVE(wka(p:q),crn(p:q,1),r)
      CALL MESSAG(cnam(31),grl+200,grt-40)
      CALL SETCLR(black) ; CALL CURVE(wka(p:q),crn(p:q,8),r)
      CALL MESSAG(cnam(32),grl+800,grt-40)
      CALL MESSAG(cnam(34),grl+400,grt+30)
      CALL ENDGRF() 
      grt=grt+300 ; grb=grb+300
      CALL NAME('Index Value','Y')   ! Axis name
      CALL tombox(y1,y2,rc,rd)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(blue)  ; CALL CURVE(wka(p:q),crn(p:q,6),r)
      CALL SETCLR(red)   ; CALL CURVE(wka(p:q),crn(p:q,2),r)
      CALL SETCLR(black) ; CALL CURVE(wka(p:q),crn(p:q,4),r)
      CALL MESSAG("b) Non-Signal-Free Chronology",grl+400,grt+30)
      CALL ENDGRF()
      grt=grt+300 ; grb=grb+300
      CALL tombox(y1,y2,rc,rd)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(blue)  ; CALL CURVE(wka(p:q),crn(p:q,7),r)
      CALL SETCLR(red)   ; CALL CURVE(wka(p:q),crn(p:q,3),r)
      CALL SETCLR(black) ; CALL CURVE(wka(p:q),crn(p:q,5),r)
      CALL MESSAG("c) Signal-Free Chronology",grl+400,grt+30) 
      CALL ENDGRF()
      grt=grt+300 ; grb=grb+300
      CALL NAME('Ratio','Y')   ! Axis name
      CALL tombox(y1,y2,re,rf)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      IF (ref1.EQ.2.OR.ref1.EQ.4.OR.ref1.EQ.5) THEN
        CALL SETCLR(cyan) ; CALL CURVE(wka(1:r),crn(1:r,50),r)
        CALL CURVE(wka(1:r),crn(1:r,52),r)
        CALL MESSAG("Expected Trend",grl+1300,grt+30)
      ENDIF 
      CALL SETCLR(blue)  ; CALL CURVE(wka(p:q),crn(p:q,36),r)
      CALL SETCLR(red)   ; CALL CURVE(wka(p:q),crn(p:q,37),r)
      CALL SETCLR(black)
      CALL MESSAG("d) Non-Signal-Free Ratios",grl+500,grt+30)
      CALL ENDGRF() 
      grt=grt+300 ; grb=grb+300
      CALL LABELS('FLOAT','X')
      CALL NAME('Calendar Year','X')     ! Axis name
      CALL tombox(y1,y2,re,rf)
      CALL SETCLR(grey) ; CALL GRID(1,1)     ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      IF (ref1.EQ.2.OR.ref1.EQ.4.OR.ref1.EQ.5) THEN
        CALL SETCLR(cyan) ; CALL CURVE(wka(1:r),crn(1:r,50),r)
        CALL CURVE(wka(1:r),crn(1:r,52),r)
        CALL MESSAG("Expected Trend",grl+1300,grt+30)
      ENDIF 
      CALL SETCLR(blue)  ; CALL CURVE(wka(p:q),crn(p:q,38),r)
      CALL SETCLR(red)   ; CALL CURVE(wka(p:q),crn(p:q,39),r)
      CALL SETCLR(black)
      CALL MESSAG("e) Signal-Free Ratios",grl+500,grt+30)
      CALL ENDGRF() ; CALL TICKS(10,'X')     ! No X ticks
      OPEN(19,FILE=jnam,IOSTAT=ios,STATUS="REPLACE")
      r=cyr(1) 
      WRITE(19,'(A20)') wnam(49)
      WRITE(19,'("Column  1 = Calendar Year")')
      WRITE(19,'("Column  2 = Count trees year")')
      IF (ref1.EQ.2.OR.ref1.EQ.4.OR.ref1.EQ.5) THEN
        WRITE(19,'("Column  3 = Mean ring +0.001")')
        WRITE(19,'("Column  4 = RCS index +0.001")')
        WRITE(19,'("Column  5 = SF RCS index +0.001")')
        WRITE(19,'("Column  6 = Mean ring ")')
        WRITE(19,'("Column  7 = RCS index ")')
        WRITE(19,'("Column  8 = SF RCS index ")')
        WRITE(19,'("Column  9 = Mean ring -0.001")')
        WRITE(19,'("Column 10 = RCS index -0.001")')
        WRITE(19,'("Column 11 = SF RCS index -0.001")')
        WRITE(19,'("Column 12 = Ratio +0.001/none")')
        WRITE(19,'("Column 13 = Ratio -0.001/none")')
        WRITE(19,'("Column 14 = SF ratio +0.001/none")')
        WRITE(19,'("Column 15 = SF ratio -0.001/none")')
        WRITE(19,'("Column 16 = Expected ratio +0.001/none")')
        WRITE(19,'("Column 17 = Expected ratio -0.001/none")')
      ELSEIF (ref1.EQ.1) THEN
        WRITE(19,'("Column  3 = Mean ring *1.0")')
        WRITE(19,'("Column  4 = RCS index *1.0")')
        WRITE(19,'("Column  5 = SF RCS index *1.0")')
        WRITE(19,'("Column  6 = Mean ring *0.7")')
        WRITE(19,'("Column  7 = RCS index *0.7")')
        WRITE(19,'("Column  8 = SF RCS index *0.7")')
        WRITE(19,'("Column  9 = Mean ring *0.4")')
        WRITE(19,'("Column 10 = RCS index *0.4")')
        WRITE(19,'("Column 11 = SF RCS index *0.4")')
        WRITE(19,'("Column 12 = Ratio 1.0/0.7")')
        WRITE(19,'("Column 13 = Ratio 0.4/0.7")')
        WRITE(19,'("Column 14 = SF ratio 1.0/0.7")')
        WRITE(19,'("Column 15 = SF ratio 0.4/0.7")')
        WRITE(19,'("Column 16 = Expected ratio 1.0/0.7")')
        WRITE(19,'("Column 17 = Expected ratio 0.4/0.7")')
      ELSEIF (ref1.NE.6) THEN
        WRITE(19,'("Column  3 = Mean ring *1.4")')
        WRITE(19,'("Column  4 = RCS index *1.4")')
        WRITE(19,'("Column  5 = SF RCS index *1.4")')
        WRITE(19,'("Column  6 = Mean ring *1.0")')
        WRITE(19,'("Column  7 = RCS index *1.0")')
        WRITE(19,'("Column  8 = SF RCS index *1.0")')
        WRITE(19,'("Column  9 = Mean ring *0.6")')
        WRITE(19,'("Column 10 = RCS index *0.6")')
        WRITE(19,'("Column 11 = SF RCS index *0.6")')
        WRITE(19,'("Column 12 = Ratio 1.4/1.0")')
        WRITE(19,'("Column 13 = Ratio 0.6/1.0")')
        WRITE(19,'("Column 14 = SF ratio 1.4/1.0")')
        WRITE(19,'("Column 15 = SF ratio 0.6/1.0")')
        WRITE(19,'("Column 16 = Expected ratio 1.4/1.0")')
        WRITE(19,'("Column 17 = Expected ratio 0.6/1.0")')
      ENDIF
      DO i=1,r
        WRITE(19,'(I4,I6,15F8.3)') i-1+cfy(1),num(i,1),crn(i,1:3), &
          crn(i,9),crn(i,6:7),crn(i,8),crn(i,4:5), &
          crn(i,36:39),crn(i,50),crn(i,52)
      ENDDO
      CLOSE(19)
      RETURN 
      END SUBROUTINE R1Fig03d
!-------------------------------------------------------------------
      SUBROUTINE R1Fig02d(ref1)  ! Plots figures
      IMPLICIT NONE                 
      INTEGER,INTENT(IN)     :: ref1 
      REAL(8),DIMENSION(mxy) :: wka
      CHARACTER(12)          :: jnam 
      REAL(8)                :: ra,rb,rc,rd  
      INTEGER                :: i=1,r,p,q
      IF     (ref1.EQ.1) THEN
        ra=0.D0 ; rb=1.49D0 ; rc=0.5D0 ; rd=1.4D0
        jnam="R1fig4.col" ; wnam(49)="Finnish trees"
      ELSEIF (ref1.EQ.2) THEN
        ra=0.D0 ; rb=1.49D0  ; rc=0.8D0 ; rd=1.19D0
        jnam="R1sm2.col" ; wnam(49)="Finnish trees"
      ELSEIF (ref1.EQ.3) THEN
        ra=0.D0 ; rb=1.4D0 ; rc=0.5D0 ; rd=1.49D0
        jnam="R1fig2.col" ; wnam(49)="Yamal trees"
      ELSEIF (ref1.EQ.4) THEN
        ra=0.D0 ; rb=1.49D0 ; rc=0.8D0 ; rd=1.2D0
        jnam="R1sm3.col" ; wnam(49)="Finnish trees"
      ELSEIF (ref1.EQ.5) THEN
        ra=0.D0 ; rb=1.7D0 ; rc=0.84D0 ; rd=1.2D0
        jnam="R1sm5.col" ; wnam(49)="Finnish trees"
      ENDIF
      p=sfy(mx) ; q=MIN(sly(mx),400) ; r=q-p+1 
      wka(1:q)=(/(DBLE(i),i=1,q)/)
      grl=200 ; grr=2000 ; grt=140 ; grb=520
      CALL TICKS(5,'X')        ! No X ticks
      CALL plot_trees(400-p+1,mcnt(p:400,12))  
      CALL LABELS('NONE','X')
      CALL NAME('','X')   ! Axis name
      CALL NAME('Ring Width','Y')   ! Axis name
      CALL tombox(1,400,ra,rb)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(blue)  ; CALL CURVE(wka(p:q),crn(p:q,16),r)
      CALL MESSAG(cnam(30),grl+1400,grt-40)
      CALL SETCLR(red)   ; CALL CURVE(wka(p:q),crn(p:q,12),r)
      CALL MESSAG(cnam(31),grl+200,grt-40)
      CALL SETCLR(black) ; CALL CURVE(wka(p:q),crn(p:q,14),r)
      CALL MESSAG(cnam(32),grl+800,grt-40)
      CALL MESSAG(cnam(33),grl+500,grt+30)
      CALL ENDGRF()
      grt=530 ; grb=910
      CALL tombox(1,400,ra,rb)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(blue)  ; CALL CURVE(wka(p:q),crn(p:q,17),r)
      CALL SETCLR(red)   ; CALL CURVE(wka(p:q),crn(p:q,13),r)
      CALL SETCLR(black) ; CALL CURVE(wka(p:q),crn(p:q,15),r)
      CALL MESSAG("b) Signal-Free RCS Curves",grl+500,grt+30) 
      CALL ENDGRF()
      grt=920 ; grb=1300
      CALL NAME('Ratio','Y')   ! Axis name
      CALL tombox(1,400,rc,rd)
      CALL SETCLR(grey) ; CALL GRID(1,1)     ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(blue)
      CALL CURVE(wka(p:q),crn(p:q,16)/crn(p:q,14),r)
      CALL SETCLR(red)
      CALL CURVE(wka(p:q),crn(p:q,12)/crn(p:q,14),r)
      CALL SETCLR(black)
      CALL MESSAG("c) Ratios Non-Signal-Free RCS",grl+500,grt+30)
      CALL ENDGRF() ; CALL TICKS(10,'X')     ! No X ticks
      grt=1310 ; grb=1690
      CALL LABELS('FLOAT','X') 
      CALL NAME('Ring Age','X')     ! Axis name
      CALL tombox(1,400,rc,rd)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(blue)
      CALL CURVE(wka(p:q),crn(p:q,17)/crn(p:q,15),r)
      CALL SETCLR(red)
      CALL CURVE(wka(p:q),crn(p:q,13)/crn(p:q,15),r)
      CALL SETCLR(black)
      CALL MESSAG("d) Ratios Signal-free RCS",grl+500,grt+30)
      CALL ENDGRF() 
      OPEN(19,FILE=jnam,IOSTAT=ios,STATUS="REPLACE")
      r=sly(mx) 
      WRITE(19,'(A20)') wnam(49)
      WRITE(19,'("Column  1 = Ring age")')
      WRITE(19,'("Column  2 = Count trees by age")')
      IF (ref1.EQ.1) THEN
        WRITE(19,'("Column  3 = Ring *1.0")')
        WRITE(19,'("Column  4 = SF ring *1.0")')
        WRITE(19,'("Column  5 = Ring *0.7")')
        WRITE(19,'("Column  6 = SF ring *0.7")')
        WRITE(19,'("Column  7 = Ring *0.4")')
        WRITE(19,'("Column  8 = SF ring *0.4")')
        WRITE(19,'("Column  9 = Ring 1.0/0.7")')
        WRITE(19,'("Column 10 = Ring 0.4/0.7")')
        WRITE(19,'("Column 11 = SF ring 1.0/0.7")')
        WRITE(19,'("Column 12 = SF ring 0.4/0.7")')
      ELSEIF (ref1.EQ.5.OR.ref1.EQ.4.OR.ref1.EQ.2) THEN
        WRITE(19,'("Column  3 = Ring +0.001")')
        WRITE(19,'("Column  4 = SF ring +0.001")')
        WRITE(19,'("Column  5 = Ring")')
        WRITE(19,'("Column  6 = SF ring")')
        WRITE(19,'("Column  7 = Ring -0.001")')
        WRITE(19,'("Column  8 = SF ring -0.001")')
        WRITE(19,'("Column  9 = Ring +0.001/none")')
        WRITE(19,'("Column 10 = Ring -0.001/none")')
        WRITE(19,'("Column 11 = SF ring +0.001/none")')
        WRITE(19,'("Column 12 = SF ring -0.001/none")')
      ELSE
        WRITE(19,'("Column  3 = Ring *1.4")')
        WRITE(19,'("Column  4 = SF ring *1.4")')
        WRITE(19,'("Column  5 = Ring *1.0")')
        WRITE(19,'("Column  6 = SF ring *1.0")')
        WRITE(19,'("Column  7 = Ring *0.6")')
        WRITE(19,'("Column  8 = SF ring *0.6")')
        WRITE(19,'("Column  9 = Ring 1.4/1.0")')
        WRITE(19,'("Column 10 = Ring 0.6/1.0")')
        WRITE(19,'("Column 11 = SF ring 1.4/1.0")')
        WRITE(19,'("Column 12 = SF ring 0.6/1.0")')
      ENDIF
      DO i=1,r
        WRITE(19,'(I4,I6,10F8.3)') i,mcnt(i,12),crn(i,12:17), &
          crn(i,12)/crn(i,14),crn(i,16)/crn(i,14), &
          crn(i,13)/crn(i,15),crn(i,17)/crn(i,15)
      ENDDO
      CLOSE(19)
      RETURN 
      END SUBROUTINE R1Fig02d
!-------------------------------------------------------------------
      SUBROUTINE R1Fig01d()  
      IMPLICIT NONE                 
      REAL(8),DIMENSION(mxy) :: wka
      INTEGER                :: i=1,p,r,q
      r=400 ; wka(1:r)=(/(DBLE(i),i=1,r)/)
      grl=200 ; grr=2000 ; grt=140 ; grb=430 
      CALL plot_trees(r-1,num(2:r,13))  
      CALL LABELS('NONE','X')
      CALL NAME('','X')     ! Axis name
      CALL NAME('Ring Width','Y')   ! Axis name
      CALL tombox(2,400,0.2D0,1.3D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(blue) ; CALL CURVE(wka(1:r),crn(1:r,13),r)
      CALL MESSAG("Mean Ring Width",grl+300,grt-45)
      CALL SETCLR(red) ; CALL CURVE(wka(1:r),crn(1:r,14),r)
      CALL MESSAG("Signal-free Ring Width",grl+1100,grt-45)
      CALL SETCLR(black)
      CALL MESSAG("a) "//wnam(3),grl+400,grt+40)
      CALL ENDGRF() 

      grt=grt+300 ; grb=grb+300
      CALL plot_trees(r-1,num(2:r,11))  
      CALL LABELS('FLOAT','X')
      CALL NAME('Ring Age','X')     ! Axis name
      CALL tombox(2,400,0.2D0,1.3D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(blue) ; CALL CURVE(wka(1:r),crn(1:r,11),r)
      CALL SETCLR(red) ; CALL CURVE(wka(1:r),crn(1:r,12),r)
      CALL SETCLR(black)
      CALL MESSAG("b) "//wnam(1),grl+400,grt+40)
      CALL ENDGRF()  

      r=cyr(3) ; wka(1:r)=(/(DBLE(i),i=cfy(3),cly(3))/)
      p=cly(3)-1200 ; grt=grt+440 ; grb=grb+440
      CALL LABELS('NONE','X')
      CALL NAME('','X')              ! Axis name
      CALL NAME('Index Value','Y')   ! Axis name
      CALL plot_trees(p,num(1:p,3)) 
      CALL tombox(1201,cly(3),0.3D0,1.6D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(blue) ; CALL CURVE(wka(1:r),crn(1:r,3),r)
      CALL MESSAG("One-curve RCS Chronology",grl+30,grt-45)
      CALL SETCLR(red)  ; CALL CURVE(wka(1:r),crn(1:r,4),r)
      CALL MESSAG("One-curve, Signal-free RCS Chronology",grl+800,grt-45)
      CALL SETCLR(black)    
      CALL MESSAG("c) "//wnam(3),grl+25,grt+30)
      CALL ENDGRF()  

      r=cyr(1) ; wka(1:r)=(/(DBLE(i),i=cfy(1),cly(1))/)
      p=1202-cfy(1) ; q=r-p+1
      grt=grt+300 ; grb=grb+300
      CALL LABELS('FLOAT','X')
      CALL NAME('Calendar Year','X')           ! Axis name
      CALL plot_trees(q,num(p:r,1))  
      CALL tombox(1201,cly(1),0.3D0,1.6D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(blue) ; CALL CURVE(wka(1:r),crn(1:r,1),r)
      CALL SETCLR(red) ; CALL CURVE(wka(1:r),crn(1:r,2),r)
      CALL SETCLR(black)    
      CALL MESSAG("d) "//Wnam(1),grl+25,grt+30)
      CALL ENDGRF()
      RETURN 
      END SUBROUTINE R1Fig01d
!-----------------------------------------------------------------------
      SUBROUTINE R1Fig01()  ! Signal free RCS subfoss & modern
      IMPLICIT NONE                 
      REAL(8) :: mn1,mn2
      INTEGER :: i,j,r,p,q    ! Tornetrask AD and Tornetrask Modern
      CALL det_default() ; idt=-2
      cnam(1)="../../raw/rcs/torn-all.raw"
      cnam(3)="../../raw/rcs/torn-liv.raw"
      wnam(1)="Tornetrask TRW"
      wnam(3)="Tornetrask Living"
      nc=0 ; CALL read_rft(cnam(1))
      cf=1 ; sfo=1 ; CALL detrend()     ! 1= Yamal not Sig-free
      r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
      num(1:r,cf+10)=mcnt(1:r,mx)       ! 11= RCS 
      cf=2 ; sfo=2 ; CALL detrend()     ! 2= Yamal Sig-free
      r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
      num(1:r,cf+10)=mcnt(1:r,mx)       ! 12= RCS
      OPEN(74,FILE="zzz.prn",IOSTAT=ios,STATUS="REPLACE")
      r=0 ; j=0 ; mn1=0.D0
      DO i=1,nc      
        IF (ly(i).GT.1950.AND.yr(i).LT.600) THEN
          p=ad(i) ; q=p+yr(i)-1 
          r=r+yr(i) ; mn1=mn1+SUM(dx(p:q))
        ENDIF
      ENDDO
      mn1=mn1/DBLE(r) 
      WRITE(74,'(2I7,F8.3)') j,r,mn1
      nc=0 ; CALL read_rft(cnam(3)) 
      cf=3 ; sfo=1 ; CALL detrend()     ! 13= Living not Sig-free
      r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
      num(1:r,cf+10)=mcnt(1:r,mx)       ! 13= RCS 
      cf=4 ; sfo=2 ; CALL detrend()     ! 14= Living Sig-free
      r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
      num(1:r,cf+10)=mcnt(1:r,mx)       ! 15= RCS
      j=0 ; r=0 ; mn2=0.D0
      DO i=1,nc      
        p=ad(i) ; q=p+yr(i)-1 
        r=r+yr(i) ; mn2=mn2+SUM(dx(p:q))
      ENDDO
      mn2=mn2/DBLE(r) 
      WRITE(74,'(2I7,F8.3)') j,r,mn2
      CLOSE(74)
      DO j=1,3,2              ! Remove count < 5
        cf=j ; r=cyr(j)
        ID: DO i=1,r ; IF (num(i,j).GE.5) EXIT ID ; ENDDO ID
        crn(1:r-i+1,j:j+1)=crn(i:r,j:j+1)
        num(1:r-i+1,j:j+1)=num(i:r,j:j+1)
        cfy(j:j+1)=cfy(j)+i-1 ; cyr(j:j+1)=cly(j)-cfy(j)+1
        CALL splinet(r,crn(1:r,j),30,crn(1:r,j)) ! Smooth chronology
        CALL splinet(r,crn(1:r,j+1),30,crn(1:r,j+1))
        IF (j.EQ.3) THEN
          crn(1:r,3)=crn(1:r,3)*mn1/mn2
          crn(1:r,4)=crn(1:r,4)*mn1/mn2
        ENDIF
      ENDDO
      r=cyr(3) ; q=cly(3)-1200  ; p=q-r+1
      num(p:q,3)=num(1:r,3) ; num(1:p-1,3)=0
      OPEN(19,FILE="R1fig1.col",IOSTAT=ios,STATUS="REPLACE")
      WRITE(19,'("Column 1 = Count living trees by age")')
      WRITE(19,'("Column 2 = Count All trees by age")')
      WRITE(19,'("Column 3 = Living mean ring width")')
      WRITE(19,'("Column 4 = Living SF ring width")')
      WRITE(19,'("Column 5 = All trees mean ring width")')
      WRITE(19,'("Column 6 = All trees SF ring width")')
      DO i=1,619
        WRITE(19,'(I4,2I6,4F8.3)') i,num(i,13),num(i,11), &
          crn(i,13),crn(i,14),crn(i,11),crn(i,12)
      ENDDO
      WRITE(19,*) ; r=cyr(4)
      WRITE(19,'("Column 1 = Count living trees by year")')
      WRITE(19,'("Column 2 = Living one-curve RCS")')
      WRITE(19,'("Column 3 = Living one-curve SF RCS")')
      DO i=1,r
        WRITE(19,'(2I6,2F8.3)') i-1+cfy(4),num(i,4),crn(i,3:4)
      ENDDO
      WRITE(19,*) ; r=cyr(1)
      WRITE(19,'("Column 1 = Count living trees by year")')
      WRITE(19,'("Column 2 = All trees one-curve RCS")')
      WRITE(19,'("Column 3 = All trees one-curve SF RCS")')
      DO i=1,r
        WRITE(19,'(2I6,2F8.3)') i-1+cfy(1),num(i,1),crn(i,1:2)
      ENDDO
      CLOSE(19)
      RETURN 
      END SUBROUTINE R1Fig01
!--------------------------------------------------------------------
      SUBROUTINE HB_rand(ref1)  ! RCS sub-fossil trees
      IMPLICIT NONE   ! u=tree age, v=chronology years               
      INTEGER,INTENT(IN)    :: ref1
      INTEGER               :: i,j,p,q,r,u,v
      REAL(8)               :: r1,r2
      INTEGER,DIMENSION(12) :: seed
      seed=10 ; CALL RANDOM_SEED(put=seed) ! Initialise random sequence
      ad(1)=1   ! Sloping line range 1.5 to 0.5
      SELECT CASE (ref1)      
      CASE (1)  ! Sub-fossil with -ve slope
        nc=81 ; u=200 ; v=1000 ; r1=+0.5D0
        wnam(20)="randb81dn.raw"  
      CASE (2)  ! Sub-fossil with +ve slope
        nc=81 ; u=200 ; v=1000 ; r1=-0.5D0
        wnam(20)="randb81up.raw"  
      CASE (3)  ! Living with -ve slope
        nc=30 ; u=100 ; v=400 ; r1=+0.5D0 
        wnam(20)="randb81dn.raw"  
      CASE (4)  ! Living with -ve slope
        nc=30 ; u=100 ; v=400 ; r1=-0.5D0   
        wnam(20)="randb81up.raw"  
      END SELECT
      DO i=1,v   ! Chronology signal
        crn(i,ref1)=1.0D0+(r1/DBLE(v))*(DBLE((v+1)/2-i))
      ENDDO
      DO i=1,nc               ! Measurement series
        p=1+(i-1)*10 ; fy(i)=p
        IF (ref1.LE.2) THEN   ! 200yr old trees
          q=p+199 
        ELSE                  ! Trees finish on last year
          q=v  
        ENDIF
        ly(i)=q ; r=q-p+1 ; yr(i)=r ; pth(i)=fy(i)-1
        WRITE(nam(i),'("Tree",I2)') i
        pthr(i)=0.1D0 ; ad(i+1)=ad(i)+r ; tre(i)=i
        DO j=1,r               ! Add 10% amplitude random noise
          CALL RANDOM_NUMBER(r1) ; r2=(9.5D0+r1)/10.0D0
          x(ad(i)-1+j)=r2*crn(p-1+j,ref1)
        ENDDO
      ENDDO 
      xok(1:ad(nc+1))=TR
!     CALL write_raw("randb30.raw",x)
!     CALL write_pith("randb30.pth")
      RETURN 
      END SUBROUTINE HB_rand
!------------------------------------------------------------------------
      SUBROUTINE R1Fig07()  ! Slope adjust no change to mean
      IMPLICIT NONE
      INTEGER :: j,r
      REAL(8) :: r1
      CALL det_default() ; idt=-2 ; cf=4 
      DO j=1,4
        CALL HB_rand(j)
        cf=cf+1 ; sfo=1 ; CALL detrend()
        r=sly(mx) ; crn(1:r,cf+20)=mval(1:r,mx)
        cyr(cf+20)=r ; num(1:r,cf+20)=mcnt(1:r,mx)
        cf=cf+1 ; sfo=2 ; CALL detrend()
        r=sly(mx) ; crn(1:r,cf+20)=mval(1:r,mx)
        num(1:r,cf+20)=mcnt(1:r,mx) ; r=cyr(cf)
        r1=SUM(crn(1:r,j)*DBLE(num(1:r,cf)))/DBLE(SUM(num(1:r,cf)))  
        crn(1:r,j)=crn(1:r,j)/r1
      ENDDO
      RETURN
      END SUBROUTINE R1Fig07
!--------------------------------------------------------------
      SUBROUTINE R1Fig07d()  ! Plots Best Fit
      IMPLICIT NONE                 
      REAL(8),DIMENSION(mxy) :: wka
      INTEGER                :: i=1,p,q,r
      p=cfy(5) ; q=cly(5) ; r=cyr(5) 
      wka(1:r)=(/(DBLE(i),i=p,q)/)   
      grl=200 ; grr=2000 ; grt=140 ; grb=490
      CALL LABELS('NONE','X')
      CALL NAME('','X')            ! Axis name
      CALL NAME('Index Value','Y') ! Axis name
      CALL plot_trees(r,num(1:r,5))  
      CALL tombox(p,q,0.7D0,1.3D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(black) ; CALL CURVE(wka(1:r),crn(1:r,1),r)
      CALL MESSAG("Expected Trend",grl+1450,grt-45)
      CALL MESSAG("a) Sub-Fossil Negative Slope",grl+600,grt+30)
      CALL SETCLR(blue)  ; CALL CURVE(wka(1:r),crn(1:r,5),r)
      CALL MESSAG("Simple RCS Chronology",grl+20,grt-45)
      CALL SETCLR(red)   ; CALL CURVE(wka(1:r),crn(1:r,6),r)
      CALL MESSAG("Signal-free RCS Chronology",grl+720,grt-45)
      CALL SETCLR(black) ; CALL ENDGRF()
      grt=grt+360 ; grb=grb+360
      CALL LABELS('FLOAT','X')
      CALL tombox(p,q,0.7D0,1.3D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(black) ; CALL CURVE(wka(1:r),crn(1:r,2),r)
      CALL MESSAG("b) Sub-Fossil Positive Slope",grl+600,grt+30) 
      CALL SETCLR(blue)  ; CALL CURVE(wka(1:r),crn(1:r,7),r)
      CALL SETCLR(red)   ; CALL CURVE(wka(1:r),crn(1:r,8),r)
      CALL SETCLR(black) ; CALL ENDGRF()
      p=cfy(9) ; q=cly(9) ; r=cyr(9) 
      wka(1:r)=(/(DBLE(i),i=p,q)/)   
      grt=grt+420 ; grb=grb+420
      CALL LABELS('NONE','X')
      CALL plot_trees(r,num(1:r,9))  
      CALL tombox(p,q,0.8D0,1.4D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(black) ; CALL CURVE(wka(1:r),crn(1:r,3),r)
      CALL MESSAG("c) Living Tree Negative Slope",grl+600,grt+30) 
      CALL SETCLR(blue)  ; CALL CURVE(wka(1:r),crn(1:r,9),r)
      CALL SETCLR(red)   ; CALL CURVE(wka(1:r),crn(1:r,10),r)
      CALL SETCLR(black) ; CALL ENDGRF()
      grt=grt+360 ; grb=grb+360
      CALL NAME('Nominal Calendar Year','X')   ! Axis name
      CALL LABELS('FLOAT','X')
      CALL tombox(p,q,0.65D0,1.2D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(black) ; CALL CURVE(wka(1:r),crn(1:r,4),r)
      CALL MESSAG("d) Living Tree Positive Slope",grl+600,grt+30) 
      CALL SETCLR(blue)  ; CALL CURVE(wka(1:r),crn(1:r,11),r)
      CALL SETCLR(red)   ; CALL CURVE(wka(1:r),crn(1:r,12),r)
      CALL SETCLR(black) ; CALL ENDGRF()
      OPEN(19,FILE="R1fig7.col",IOSTAT=ios,STATUS="REPLACE")
      r=1000 
      WRITE(19,'("Artificial trees")') 
      WRITE(19,'("Column   1 = Nominal calendar year")')
      WRITE(19,'("Column   2 = Count sub-fossil")')
      WRITE(19,'("Column   3 = Count living")')
      WRITE(19,'("Column   4 = Sub-fossil negative expected")')
      WRITE(19,'("Column   5 = Sub-fossil positive expected")')
      WRITE(19,'("Column   6 = Living negative expected")')
      WRITE(19,'("Column   7 = living positive expected")')
      WRITE(19,'("Column   8 = Sub-fossil negative simple RCS")')
      WRITE(19,'("Column   9 = Sub-fossil negative SF RCS")')
      WRITE(19,'("Column  10 = Sub-fossil positive simple RCS")')
      WRITE(19,'("Column  11 = Sub-fossil positive SF RCS")')
      WRITE(19,'("Column  12 = Living negative simple RCS")')
      WRITE(19,'("Column  13 = Living negative SF RCS")')
      WRITE(19,'("Column  14 = Living positive simple RCS")')
      WRITE(19,'("Column  15 = Living positive SF RCS")')
      DO i=1,r
        WRITE(19,'(I4,2I6,12F8.3)') i,num(i,5),num(i,9),crn(i,1:12)
      ENDDO
      CLOSE(19)
      RETURN 
      END SUBROUTINE R1Fig07d
!-------------------------------------------------------------------
      SUBROUTINE R1Fig08d()  ! Plots Best Fit
      IMPLICIT NONE                 
      REAL(8),DIMENSION(mxy) :: wka
      INTEGER                :: i=1,r
      r=cyr(25) ; wka(1:500)=(/(DBLE(i),i=1,500)/)   
      grl=200 ; grr=2000 ; grt=140 ; grb=490
      CALL LABELS('NONE','X')
      CALL NAME('','X')            ! Axis name
      CALL NAME('Ring Width','Y') ! Axis name
      CALL plot_trees(r,num(1:r,25))  
      CALL tombox(1,r,0.7D0,1.3D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(black)
      CALL MESSAG("a) Sub-Fossil Negative Slope",grl+600,grt+30)
      CALL SETCLR(blue)  ; CALL CURVE(wka(2:r),crn(2:r,25),r-1)
      CALL MESSAG("Simple RCS Curve",grl+100,grt-45)
      CALL SETCLR(red)   ; CALL CURVE(wka(2:r),crn(2:r,26),r-1)
      CALL MESSAG("Signal-free RCS Curve",grl+800,grt-45)
      CALL SETCLR(black) ; CALL ENDGRF()
      r=cyr(27) ; grt=grt+360 ; grb=grb+360
      CALL LABELS('FLOAT','X')
      CALL tombox(1,r,0.7D0,1.3D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(black) 
      CALL MESSAG("b) Sub-Fossil Positive Slope",grl+600,grt+30) 
      CALL SETCLR(blue)  ; CALL CURVE(wka(2:r),crn(2:r,27),r-1)
      CALL SETCLR(red)   ; CALL CURVE(wka(2:r),crn(2:r,28),r-1)
      CALL SETCLR(black) ; CALL ENDGRF()
      r=cyr(29) ; grt=grt+420 ; grb=grb+420
      CALL LABELS('NONE','X')
      CALL plot_trees(r,num(1:r,29))  
      CALL tombox(1,r,0.7D0,1.3D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(black) 
      CALL MESSAG("c) Living Tree Negative Slope",grl+600,grt+30) 
      CALL SETCLR(blue)  ; CALL CURVE(wka(2:r),crn(2:r,29),r-1)
      CALL SETCLR(red)   ; CALL CURVE(wka(2:r),crn(2:r,30),r-1)
      CALL SETCLR(black) ; CALL ENDGRF()
      r=cyr(31) ; grt=grt+360 ; grb=grb+360
      CALL NAME('Ring Age','X')   ! Axis name
      CALL LABELS('FLOAT','X')
      CALL tombox(1,r+1,0.7D0,1.3D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(black) 
      CALL MESSAG("d) Living Tree Positive Slope",grl+600,grt+30) 
      CALL SETCLR(blue)  ; CALL CURVE(wka(2:r),crn(2:r,31),r-1)
      CALL SETCLR(red)   ; CALL CURVE(wka(2:r),crn(2:r,32),r-1)
      CALL SETCLR(black) ; CALL ENDGRF()
      OPEN(19,FILE="R1fig6.col",IOSTAT=ios,STATUS="REPLACE")
      r=410 
      WRITE(19,'("Artificial trees")') 
      WRITE(19,'("Column   1 = Ring age")')
      WRITE(19,'("Column   2 = Count sub-fossil")')
      WRITE(19,'("Column   3 = Count living")')
      WRITE(19,'("Column   4 = Sub-fossil negative simple RCS")')
      WRITE(19,'("Column   5 = Sub-fossil negative SF RCS")')
      WRITE(19,'("Column   6 = Sub-fossil positive simple RCS")')
      WRITE(19,'("Column   7 = Sub-fossil positive SF RCS")')
      WRITE(19,'("Column   8 = Living negative simple RCS")')
      WRITE(19,'("Column   9 = Living negative SF RCS")')
      WRITE(19,'("Column  10 = Living positive simple RCS")')
      WRITE(19,'("Column  11 = Living positive SF RCS")')
      DO i=1,r
        WRITE(19,'(I4,2I6,8F8.3)') i,num(i,25),num(i,29),crn(i,25:32)
      ENDDO
      CLOSE(19)
      RETURN 
      END SUBROUTINE R1Fig08d
!-------------------------------------------------------------------
      SUBROUTINE RC_sm21(ref1)  ! Reduced rings final years 
      IMPLICIT NONE
      INTEGER,INTENT(IN)   :: ref1
      INTEGER              :: i,j,p,q,r
      CHARACTER(12)        :: jnam
      REAL(8),DIMENSION(3) :: fac
      OPEN(74,FILE="zzz.prn",IOSTAT=ios,STATUS="REPLACE")
      CALL det_default() ; idt=-2
      fac=(/1.D0,0.7D0,0.4D0/) ; cf=1 
      IF (ref1.EQ.1) THEN
        cnam(18)="../../raw/rcs/finn100.raw"
!       jnam="R1sm7-11.col"
      ELSEIF (ref1.EQ.2) THEN
        cnam(18)="../../raw/rcs/finn150.raw"
        jnam="R1sm7.col"
      ELSEIF (ref1.EQ.3) THEN
        cnam(18)="../../raw/rcs/finn200.raw"
        jnam="R1sm8.col"
      ELSEIF (ref1.EQ.4) THEN
        cnam(18)="../../raw/rcs/finn250.raw"
        jnam="R1sm9.col"
      ELSEIF (ref1.EQ.5) THEN
        cnam(18)="../../raw/rcs/finn300.raw"
        jnam="R1sm10.col"
      ELSEIF (ref1.EQ.6) THEN
        cnam(18)="../../raw/rcs/finn350.raw"
        jnam="R1sm11.col"
      ENDIF
      DO j=1,3
        nc=0 ; CALL read_rft(cnam(18))
        DO i=1,nc
          p=MAX(fy(i),1920) ; q=ly(i)
          IF (p.LT.q) THEN
            p=ad(i)+(p-fy(i)) ; q=ad(i+1)-1
            x(p:q)=x(p:q)*fac(j)
          ENDIF
        ENDDO
        IF (j.EQ.2) THEN
          cf=cf+5 ; idt=0 ; CALL detrend() ; cf=cf-5  ! cf=8
        ELSEIF (j.EQ.3) THEN
          cf=cf+4 ; idt=0 ; CALL detrend() ; cf=cf-4  ! cf=9
        ELSE
          idt=0 ; CALL detrend()         ! cf=1
        ENDIF
        idt=-2 ; cf=cf+1 ; sfo=1 ; CALL detrend()
        r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
        num(1:r,cf+10)=mcnt(1:r,mx)
        cf=cf+1 ; sfo=2 ; CALL detrend()
        r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
        num(1:r,cf+10)=mcnt(1:r,mx)
        WRITE(74,'(4I6)') j,cfy(cf-2),cly(cf-2),cyr(cf-2)
        WRITE(74,'(4I6)') cf,cfy(cf-1),cly(cf-1),cyr(cf-1)
        WRITE(74,'(4I6)') cf,cfy(cf),cly(cf),cyr(cf)
      ENDDO
      r=cyr(1)
      crn(1:r,36)=crn(1:r,6)/crn(1:r,4)
      crn(1:r,37)=crn(1:r,2)/crn(1:r,4)
      crn(1:r,38)=crn(1:r,7)/crn(1:r,5)
      crn(1:r,39)=crn(1:r,3)/crn(1:r,5)
      DO i=1,12
        WRITE(74,'(4I6)') i,cfy(i),cly(i),cyr(i)
      ENDDO
      CLOSE(74)
      cnam(30)="1920+ Rings *0.4"
      cnam(31)="Measured Values"
      cnam(32)="1920+ Rings *0.7"
      cnam(33)="a) Finnish Trees - RCS Curves"
      cnam(34)="a) Finnish Trees - Ring-Width"

      OPEN(19,FILE=jnam,IOSTAT=ios,STATUS="REPLACE")
      r=cyr(1) 
      WRITE(19,'(A20)') cnam(18)(15:25)//" trees"
      WRITE(19,'("Column  1 = Calendar Year")')
      WRITE(19,'("Column  2 = Count trees year")')
      WRITE(19,'("Column  3 = Mean ring *1.0")')
      WRITE(19,'("Column  4 = RCS index *1.0")')
      WRITE(19,'("Column  5 = SF RCS index *1.0")')
      WRITE(19,'("Column  6 = Mean ring *0.7")')
      WRITE(19,'("Column  7 = RCS index *0.7")')
      WRITE(19,'("Column  8 = SF RCS index *0.7")')
      WRITE(19,'("Column  9 = Mean ring *0.4")')
      WRITE(19,'("Column 10 = RCS index *0.4")')
      WRITE(19,'("Column 11 = SF RCS index *0.4")')
      WRITE(19,'("Column 12 = Ratio 1.0/0.7")')
      WRITE(19,'("Column 13 = Ratio 0.4/0.7")')
      WRITE(19,'("Column 14 = SF ratio 1.0/0.7")')
      WRITE(19,'("Column 15 = SF ratio 0.4/0.7")')
      DO i=1,r
        WRITE(19,'(I4,I6,15F8.3)') i-1+cfy(1),num(i,1),crn(i,1:3), &
          crn(i,9),crn(i,6:7),crn(i,8),crn(i,4:5),crn(i,36:39)
      ENDDO
      CLOSE(19)
      RETURN
      END SUBROUTINE RC_sm21
!--------------------------------------------------------------
      SUBROUTINE R1SM01d()  
      IMPLICIT NONE                 
      REAL(8),DIMENSION(mxy)    :: wka
      INTEGER                   :: i=1,r
      INTEGER,DIMENSION(7)      :: col
      CHARACTER(3),DIMENSION(7) :: lab
      col=(/blue,cyan,red,lime,yellow,brown,black/)
      lab=(/"0.4","0.5","0.6","0.7","0.8","0.9","1.0"/)
      r=cyr(1) ; wka(1:r)=(/(DBLE(i),i=cfy(1),cly(1))/)   
      grl=200 ; grr=2000 ; grt=140 ; grb=430
      CALL TICKS(0,'X')            ! No X ticks
      CALL LABELS('NONE','X') 
      CALL NAME('','X')     ! Axis name
      CALL plot_trees(r,num(1:r,1))  
      CALL NAME('Index Value','Y')   ! Axis name
      CALL tombox(cfy(1),cly(1),0.5D0,1.4D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      DO i=1,7
        CALL SETCLR(col(i))
        CALL CURVE(wka(1:r),crn(1:r,i),r)
        CALL MESSAG(lab(i),grl+140+i*200,grt-45)
      ENDDO
      CALL SETCLR(black)
      CALL MESSAG("Step Factor",grl+10,grt-45)
      CALL MESSAG("a) Non-Signal-Free Chronology",grl+50,grt+18)
      CALL ENDGRF()
      grt=grt+300 ; grb=grb+300
      CALL tombox(cfy(1),cly(1),0.5D0,1.4D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      DO i=1,7
        CALL SETCLR(col(i))
        CALL CURVE(wka(1:r),crn(1:r,i+7),r)
      ENDDO
      CALL SETCLR(black)
      CALL MESSAG("b) Signal-Free Chronology",grl+50,grt+18) 
      CALL ENDGRF()
      grt=grt+300 ; grb=grb+300
      CALL NAME('Ratio','Y')   ! Axis name
      CALL tombox(cfy(1),cly(1),0.6D0,1.32D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      DO i=1,7
        CALL SETCLR(col(i))
        CALL CURVE(wka(1:r),crn(1:r,i+21),r)
      ENDDO
      CALL SETCLR(black)
      CALL MESSAG("c) Simple Ratio to 0.7",grl+50,grt+18) 
      CALL ENDGRF()
      grt=grt+300 ; grb=grb+300
      CALL LABELS('FLOAT','X') 
      CALL NAME('Years','X')     ! Axis name
      CALL tombox(cfy(1),cly(1),0.6D0,1.32D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      DO i=1,7
        CALL SETCLR(col(i))
        CALL CURVE(wka(1:r),crn(1:r,i+28),r)
      ENDDO
      CALL SETCLR(black)
      CALL MESSAG("d) Signal-Free Ratio to 0.7",grl+50,grt+18) 
      CALL ENDGRF()
      CALL TICKS(5,'X')            ! No X ticks
      RETURN 
      END SUBROUTINE R1SM01d
!-------------------------------------------------------------------
      SUBROUTINE R1SM01()  ! Reduced rings final years 
      IMPLICIT NONE
      INTEGER  :: i,j,p,q,r
      CALL det_default() ; idt=-2 ; cf=0 ! ; idb=2
      DO j=4,10
        nc=0 ; CALL read_rft("../../raw/rcs/finnmrg.raw")
        DO i=1,nc
          p=MAX(fy(i),1920) ; q=ly(i)
          IF (p.LT.q) THEN
            p=ad(i)+(p-fy(i)) ; q=ad(i+1)-1
            x(p:q)=x(p:q)*DBLE(j)/10.D0
          ENDIF
        ENDDO
        cf=j-3 ; sfo=1 ; CALL detrend()
        cf=j+4 ; sfo=2 ; CALL detrend()
      ENDDO
      r=cyr(1) ; p=1920-cfy(1) 
      DO j=1,7
        crn(1:r,j+21)=crn(1:r,j)/crn(1:r,4)
        crn(1:r,j+28)=crn(1:r,j+7)/crn(1:r,11)
        crn(1:r,j+14)=crn(1:r,j)/crn(1:r,j+7)
      ENDDO
      DO j=1,7
        CALL splinet(r,crn(1:r,j),50,crn(1:r,j))    
        CALL splinet(r,crn(1:r,j+7),50,crn(1:r,j+7))    
      ENDDO
      OPEN(19,FILE="R1sm1.col",STATUS="REPLACE")
      WRITE(19,'(A20)') "Finnish trees"
      WRITE(19,'("Column  1 = Calendar year")')
      WRITE(19,'("Column  2 = Count of trees")')
      WRITE(19,'("Column  3 = Simple chronology 0.4")')
      WRITE(19,'("Column  4 = Simple chronology 0.5")')
      WRITE(19,'("Column  5 = Simple chronology 0.6")')
      WRITE(19,'("Column  6 = Simple chronology 0.7")')
      WRITE(19,'("Column  7 = Simple chronology 0.8")')
      WRITE(19,'("Column  8 = Simple chronology 0.9")')
      WRITE(19,'("Column  9 = Simple chronology 1.0")')
      WRITE(19,'("Column 10 = SF chronology 0.4")')
      WRITE(19,'("Column 11 = SF chronology 0.5")')
      WRITE(19,'("Column 12 = SF chronology 0.6")')
      WRITE(19,'("Column 13 = SF chronology 0.7")')
      WRITE(19,'("Column 14 = SF chronology 0.8")')
      WRITE(19,'("Column 15 = SF chronology 0.9")')
      WRITE(19,'("Column 16 = SF chronology 1.0")')
      WRITE(19,'("Column 17 = Simple ratio 0.4")')
      WRITE(19,'("Column 18 = Simple ratio 0.5")')
      WRITE(19,'("Column 19 = Simple ratio 0.6")')
      WRITE(19,'("Column 20 = Simple ratio 0.7")')
      WRITE(19,'("Column 21 = Simple ratio 0.8")')
      WRITE(19,'("Column 22 = Simple ratio 0.9")')
      WRITE(19,'("Column 23 = Simple ratio 1.0")')
      WRITE(19,'("Column 24 = SF ratio 0.4")')
      WRITE(19,'("Column 25 = SF ratio 0.5")')
      WRITE(19,'("Column 26 = SF ratio 0.6")')
      WRITE(19,'("Column 27 = SF ratio 0.7")')
      WRITE(19,'("Column 28 = SF ratio 0.8")')
      WRITE(19,'("Column 29 = SF ratio 0.9")')
      WRITE(19,'("Column 30 = SF ratio 1.0")')
      DO i=1,r
        WRITE(19,'(2I4,28F7.3)') i-1+cfy(1),num(i,1), &
          crn(i,1:14),crn(i,22:35)
      ENDDO
      CLOSE(19)
      RETURN
      END SUBROUTINE R1SM01
!--------------------------------------------------------------
      SUBROUTINE R1Fig09(ref1)  ! Reduced rings final years 
      IMPLICIT NONE
      INTEGER,INTENT(IN)   :: ref1
      INTEGER              :: i,j,k,p,q,r,u,v
      REAL(8),DIMENSION(3) :: fac
      CALL det_default() ; idt=-2 ; cf=0
      IF (ref1.EQ.1) THEN
        fac=(/1.D0,0.7D0,0.4D0/) 
        cnam(26)="Step Decrease/No Step"
        cnam(27)="Step Increase/No Step"
      ELSE
        fac=(/+0.001D0,0.D0,-0.001D0/)
        cnam(26)="Slope Decrease/Original"
        cnam(27)="Slope Increase/Original"
      ENDIF
      cnam(1)="../../raw/rcs/finn150.raw"
      cnam(2)="../../raw/rcs/finn200.raw"
      cnam(3)="../../raw/rcs/finn250.raw"
      cnam(4)="../../raw/rcs/finn300.raw"
      cnam(5)="../../raw/rcs/finn350.raw"
      cnam(21)="a) Trees Start 1704 to 1946"
      cnam(22)="b) Trees Start 1704 to 1886"
      cnam(23)="c) Trees Start 1704 to 1774"
      cnam(24)="d) Trees Start 1620 to 1774"
      cnam(25)="e) Trees Start 1550 to 1774"
      OPEN(74,FILE="zzz.prn",IOSTAT=ios,STATUS="REPLACE")
      WRITE(74,'(" Trees  From    To   min   max range")')
      DO k=1,5
        DO j=1,3
          cf=cf+1 ; nc=0 ; CALL read_rft(cnam(k))
          IF (ref1.EQ.1) THEN
            DO i=1,nc
              p=MAX(fy(i),1920) ; q=ly(i)
              IF (p.LT.q) THEN
                p=ad(i)+(p-fy(i)) ; q=ad(i+1)-1
                x(p:q)=x(p:q)*fac(j)
              ENDIF
            ENDDO
            CALL detrend()
          ELSE
            CALL det_crnfy() ; r=cyr(cf)
            crn(1:r,49+cf)=(/(1.D0+DBLE(i)*fac(j),i=-r/2,r-r/2-1)/)
            DO i=1,nc
              r=yr(i) ; p=ad(i) ; q=p+r-1 
              u=fy(i)-cfy(cf)+1 ; v=u+r-1
              x(p:q)=x(p:q)*crn(u:v,49+cf)
            ENDDO
            CALL detrend()
          ENDIF
        ENDDO
        WRITE(74,'(6I6)') nc,MINVAL(fy(1:nc)),MAXVAL(fy(1:nc)), &
          MINVAL(yr(1:nc)),MAXVAL(yr(1:nc)),MAXVAL(yr(1:nc))-MINVAL(yr(1:nc))+1
      ENDDO
      DO i=1,15,3
        r=cyr(i) ; p=cfy(i)-cfy(13)+1 ; q=p+r-1 
        num(1:500,i+15)=0 ; num(p:q,i+15)=num(1:r,i) 
        WHERE (crn(1:r,i+1).GT.0.0001D0)
          crn(p:q,i+15)=crn(1:r,i)/crn(1:r,i+1)
          crn(p:q,i+16)=crn(1:r,i+2)/crn(1:r,i+1)
        ELSEWHERE
          crn(p:q,i+15)=0.D0
        END WHERE
        IF (ref1.EQ.2) THEN
          crn(p:q,i+49)=crn(1:r,i+49) 
          crn(p:q,i+51)=crn(1:r,i+51) 
          J1: DO j=p,q ; IF (crn(j,i+15).GT.1.D0) EXIT J1 ; ENDDO J1
          crn(p:q,i+49)=crn(p:q,i+49)/crn(j,i+49)
          J2: DO j=p,q ; IF (crn(j,i+16).LT.1.D0) EXIT J2 ; ENDDO J2
          crn(p:q,i+51)=crn(p:q,i+51)/crn(j,i+51)
        ENDIF
      ENDDO
      CLOSE(74)
      RETURN
      END SUBROUTINE R1Fig09
!--------------------------------------------------------------
      SUBROUTINE R1Fig09d(ref1)  
      IMPLICIT NONE                 
      INTEGER,INTENT(IN)     :: ref1
      REAL(8),DIMENSION(mxy) :: wka
      INTEGER                :: i=1,j,k,p,q,r,w
      r=cyr(13) ; wka(1:r)=(/(DBLE(i),i=cfy(13),cly(13))/)   
      grl=200 ; grr=2000 ; grt=140 ; grb=440
      CALL TICKS(0,'X')          ! No X ticks
      CALL LABELS('NONE','X')
      CALL NAME('','X')          ! Axis name
      CALL NAME('Ratio','Y')     ! Axis name
      DO j=1,5
        IF (j.EQ.5) THEN
          CALL LABELS('FLOAT','X')
          CALL NAME('Calendar Year','X')     ! Axis name
        ENDIF 
        i=j*3+13 ; CALL plot_treesq(r,num(1:r,i),99)  
        IF (ref1.EQ.1) THEN
          CALL tombox(cfy(13),cly(13),0.6D0,1.3D0)
        ELSE
          CALL tombox(cfy(13),cly(13),0.7D0,1.39D0)
        ENDIF
        CALL SETCLR(grey) ; CALL GRID(1,1)     ! GRIDLINES
        CALL LINWID(1) ; CALL HEIGHT(22)
        p=cfy(j*3)-cfy(13)+1 ; q=cly(j*3)-cfy(13)+1 ; w=q-p+1
        IF (j.GE.4) THEN ; p=p+3 ; w=w-3 ; ENDIF
        IF (ref1.EQ.2) THEN
          CALL SETCLR(cyan) ; CALL CURVE(wka(p:q),crn(p:q,34+i),w)
          CALL CURVE(wka(p:q),crn(p:q,36+i),w)
          IF (j.EQ.1) CALL MESSAG("Expected Trend",grl+1450,grt-40)
        ENDIF
        CALL SETCLR(blue) ! ; CALL CURVE(wka(p:q),crn(p:q,i),w)
        DO k=p,q-2,5 ; CALL CURVE(wka(k:k+2),crn(k:k+2,i+1),3) ; ENDDO
        IF (j.EQ.1) CALL MESSAG(cnam(26),grl+750,grt-40)
        CALL SETCLR(red)   ; CALL CURVE(wka(p:q),crn(p:q,i),w)
        IF (j.EQ.1) CALL MESSAG(cnam(27),grl+50,grt-40)
        CALL SETCLR(black) ; CALL MESSAG(cnam(20+j),grl+20,grt+20)
        CALL ENDGRF() ; 
        grt=grt+310 ; grb=grb+310
      ENDDO
      CALL TICKS(10,'X')     ! No X ticks
      IF (ref1.EQ.1) THEN
        OPEN(19,FILE="R1fig9.col",IOSTAT=ios,STATUS="REPLACE")
      ELSE
        OPEN(19,FILE="R1sm12.col",IOSTAT=ios,STATUS="REPLACE")
      ENDIF
      r=cyr(13) 
      DO i=21,25 ; WRITE(19,'(A)') cnam(i) ; ENDDO
      WRITE(19,'("Artificial trees")') 
      WRITE(19,'("Column   1 = Calendar year")')
      WRITE(19,'("Column   2 = Count a)")')
      WRITE(19,'("Column   3 = Count b)")')
      WRITE(19,'("Column   4 = Count c)")')
      WRITE(19,'("Column   5 = Count d)")')
      WRITE(19,'("Column   6 = Count e)")')
      WRITE(19,'("Column   7 = Increase/none a)")')
      WRITE(19,'("Column   8 = Decrease/none a)")')
      WRITE(19,'("Column   9 = Increase/none b)")')
      WRITE(19,'("Column  10 = Decrease/none b)")')
      WRITE(19,'("Column  11 = Increase/none c)")')
      WRITE(19,'("Column  12 = Decrease/none c)")')
      WRITE(19,'("Column  13 = Increase/none d)")')
      WRITE(19,'("Column  14 = Decrease/none d)")')
      WRITE(19,'("Column  15 = Increase/none e)")')
      WRITE(19,'("Column  16 = Decrease/none e)")')
      IF (ref1.EQ.2) THEN
        WRITE(19,'("Column  17 = Expected increase/none a)")')
        WRITE(19,'("Column  18 = Expected decrease/none a)")')
        WRITE(19,'("Column  19 = Expected increase/none b)")')
        WRITE(19,'("Column  20 = Expected decrease/none b)")')
        WRITE(19,'("Column  21 = Expected increase/none c)")')
        WRITE(19,'("Column  22 = Expected decrease/none c)")')
        WRITE(19,'("Column  23 = Expected increase/none d)")')
        WRITE(19,'("Column  24 = Expected decrease/none d)")')
        WRITE(19,'("Column  25 = Expected increase/none e)")')
        WRITE(19,'("Column  26 = Expected decrease/none e)")')
        DO i=1,r
          WRITE(19,'(6I4,20F7.3)') i,num(i,16),num(i,19),num(i,22), &
            num(i,25),num(i,28),crn(i,16:17),crn(i,19:20), &
            crn(i,22:23),crn(i,25:26),crn(i,28:29),crn(i,50), &
            crn(i,52:53),crn(i,55:56),crn(i,58:59), &
            crn(i,61:62),crn(i,64)
        ENDDO
      ELSE 
        DO i=1,r
          WRITE(19,'(6I4,10F7.3)') i,num(i,16),num(i,19),num(i,22), &
            num(i,25),num(i,28),crn(i,16:17),crn(i,19:20), &
            crn(i,22:23),crn(i,25:26),crn(i,28:29)
        ENDDO
      ENDIF 
      CLOSE(19)
      RETURN 
      END SUBROUTINE R1Fig09d
!-------------------------------------------------------------------
      SUBROUTINE R2Fig3d()  
      IMPLICIT NONE                 
      REAL(8),DIMENSION(mxy) :: wka
      INTEGER                :: i=1,r
      r=1001 ; wka(1:r)=(/(DBLE(i)/100.D0,i=-500,500)/) 
      grl=200 ; grr=2200 ; grt=200 ; grb=1000
      CALL LABDIG(-1,'Y')          ! Integer labels
      CALL LABDIG(1,'X')           ! 1 dec place
      CALL TICKS(2,'Y')            
      CALL NAME('z-score','X')     ! Axis name
      CALL NAME('Count of Values','Y')   ! Axis name
      CALL tombox(-3,+3,1.D0,520.D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(red) ; CALL CURVE(wka(1:r),crn(1:r,13),r)
      CALL MESSAG("RCS Indices",grl+50,grt-40)
      CALL SETCLR(cyan) ; CALL CURVE(wka(1:r),crn(1:r,14),r)
      CALL MESSAG("RCS PTrans",grl+550,grt-40)
      CALL SETCLR(blue) ; CALL CURVE(wka(1:r),crn(1:r,12),r)
      CALL MESSAG("30yr Spline",grl+1050,grt-40)
      CALL SETCLR(black)
      DO i=1,r-2,5 ; CALL CURVE(wka(i:i+2),crn(i:i+2,15),3) ; ENDDO
      CALL MESSAG("RCS Normal",grl+1550,grt-40)
      CALL LABDIG(-2,'Y')           ! Dislin labels
      CALL LABDIG(-1,'X')           ! Integer labels
      CALL TICKS(5,'Y')             
      CALL ENDGRF()
      RETURN 
      END SUBROUTINE R2Fig3d
!-------------------------------------------------------------------
      SUBROUTINE R2Fig3()  ! Distribution of various indices 
      IMPLICIT NONE
      INTEGER  :: i,j,k,m
      REAL(8)  :: rr,mn,sd
      cnam(1)="../../raw/rcs/yml-all.raw"
      nc=0 ; cf=1 ; CALL read_rft(cnam(1))
      CALL det_default()
      m=ad(nc)+yr(nc)-1 ; rr=DBLE(m)     
      crn(1:1000,10:15)=0.D0 
      DO i=1,4
        cf=i
        IF (i.EQ.1) THEN
          idt=30                 ! Spline detrend
          CALL detrend()             
        ELSEIF (i.EQ.2) THEN
          idt=-2                 ! RCS detrend
          CALL detrend()             
        ELSEIF (i.EQ.3) THEN
          itn=2 ; ind=2 ; sfo=1  ! RCS with power transform
          CALL detrend()             
        ELSEIF (i.EQ.4) THEN
          CALL randnorm(m,dx)    ! Random number sequence 
        ENDIF 
        mn=SUM(dx(1:m))/rr       ! Normalise indices
        sd=SQRT(SUM((dx(1:m)-mn)**2)/(rr-1.D0))
        dx(1:m)=(dx(1:m)-mn)/sd
        DO k=1,m                 ! Accumulate distribution
          j=NINT(dx(k)*100.D0)+501 ; j=MAX(1,MIN(1000,j))
          crn(j,11+i)=crn(j,11+i)+1.D0
        ENDDO                    ! Smooth distribution curve
        CALL splinet(1000,crn(1:1000,11+i),60,crn(1:1000,11+i))    
      ENDDO
      OPEN(19,FILE="R2fig3.col",IOSTAT=ios,STATUS="REPLACE")
      WRITE(19,'("Column  1 = Year")')
      WRITE(19,'("Column  2 = 30yr Spline")')
      WRITE(19,'("Column  3 = RCS Indices")')
      WRITE(19,'("Column  4 = RCS PTrans")')
      WRITE(19,'("Column  5 = RCS Normal")')
      DO i=1,1000
        WRITE(19,'(F8.3,4F8.2)') DBLE(i-501)/100.D0,crn(i,12:15)
      ENDDO
      CLOSE(19)
      RETURN
      END SUBROUTINE R2Fig3
!--------------------------------------------------------------
      SUBROUTINE R2Fig2d()  ! Distribution of climate 
      IMPLICIT NONE                 
      CHARACTER(32),DIMENSION(3:5) :: lab
      REAL(8),DIMENSION(mxy)       :: wka
      INTEGER :: i=1,j,r
      INTEGER :: ra
      lab=(/"b) S.Dev Sorted by Calendar Year", &
            "c) S.Dev Sorted on Index        ", &
            "d) S.Dev Scaled by Index        "/)
      r=cyr(1) ; wka(1:r)=(/(DBLE(i),i=cfy(1),cly(1))/) 
      grl=200 ; grr=2000 ; ra=140 ; grt=ra ; grb=ra+250
      CALL plot_trees(r,num(1:r,2))  
      CALL LABELS('NONE','X')
      CALL TICKS(2,'Y')            
      CALL TICKS(5,'X')            
      CALL NAME('','X')        ! Axis name
      CALL NAME('Index','Y')   ! Axis name
      CALL tombox(cfy(1),cly(1),0.1D0,2.49D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(black)
      CALL CURVE(wka(1:r),crn(1:r,2),r)
      CALL MESSAG("a) One-curve RCS Chronology",grl+50,grt+27)
      CALL ENDGRF()  
      DO j=3,5
        ra=ra+260 ; grt=ra ; grb=ra+240
        CALL tombox(cfy(1),cly(1),0.D0,1.45D0)
        CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(silver) 
        DO i=1,r
          CALL RLINE(wka(i),0.D0,wka(i),crn(i,j))
        ENDDO
        CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
        CALL SETCLR(black) ; CALL MESSAG(lab(j),grl+50,grt+30)
        CALL ENDGRF() 
      ENDDO
      ra=ra+260 ; grt=ra ; grb=ra+250
      CALL NAME('z-scores','Y')   ! Axis name
      CALL plot_trees(r,num(1:r,2))  
      CALL tombox(cfy(1),cly(1),-1.99D0,1.99D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(black)
      CALL CURVE(wka(1:r),crn(1:r,9),r)
      CALL MESSAG("e) One-curve, Normal RCS Chronology",grl+50,grt+27)
      CALL ENDGRF()
      ra=ra+260 ; grt=ra ; grb=ra+240
      CALL tombox(cfy(1),cly(1),0.D0,1.45D0)
      CALL SETCLR(silver)
      DO i=1,r
        CALL RLINE(wka(i),0.001D0,wka(i),crn(i,10))
      ENDDO
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(black) 
      CALL MESSAG("f) S.Dev Sorted by Calendar Year",grl+50,grt+30)
      CALL ENDGRF() 
      ra=ra+260 ; grt=ra ; grb=ra+240
      CALL LABELS('FLOAT','X') 
      CALL NAME('Years','X')     ! Axis name
      CALL tombox(cfy(1),cly(1),0.D0,1.45D0)
      CALL SETCLR(silver) 
      DO i=1,r
        CALL RLINE(wka(i),0.001D0,wka(i),crn(i,11))
      ENDDO
      CALL SETCLR(grey) ; CALL GRID(1,1)  ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(black) 
      CALL MESSAG("g) S.Dev Sorted on Index",grl+50,grt+30)
      CALL LABELS('FLOAT','Y') 
      CALL TICKS(5,'Y')            
      CALL TICKS(10,'X')            
      CALL ENDGRF()  
      RETURN 
      END SUBROUTINE R2Fig2d
!-------------------------------------------------------------------
      SUBROUTINE R2Fig2()  ! BFM and distribution adjust
      IMPLICIT NONE
      INTEGER                :: i,j,m,p,q,r,u,v
      INTEGER,DIMENSION(mxd) :: nx,ny   
      nc=0 ; CALL read_rft("../../raw/rcs/yml-all.raw")
      CALL det_default()
      idt=-2 ; cf=1 ; CALL detrend() ; r=cyr(1)
      crn(1:r,2:8)=0.D0 ; num(1:r,2)=0
      DO i=1,nc              ! Simple RCS
        p=ad(i) ; j=yr(i) ; q=p+j-1 ; u=fy(i)-cfy(1)+1 ; v=u+j-1
        WHERE (xok(p:q))
          num(u:v,2)=num(u:v,2)+1
          crn(u:v,3)=crn(u:v,3)+dx(p:q)     ! Sum indices
          crn(u:v,4)=crn(u:v,4)+dx(p:q)**2  ! Sum squares
        END WHERE
      ENDDO                  ! Ignore counts < 4
      JD: DO j=1,r ; IF (num(j,2).GE.4) EXIT JD ; ENDDO JD
      crn(1:r-j+1,1:4)=crn(j:r,1:4)
      num(1:r-j+1,2)=num(j:r,2)  
      cfy(1)=cfy(1)+j-1 ; r=r-j+1 ; cyr(1)=r
      crn(1:r,2)=crn(1:r,3)/DBLE(num(1:r,2))      ! CRN
      crn(1:r,3)=SQRT(MAX(crn(1:r,4)-crn(1:r,2)* &  
        crn(1:r,3),0.001D0)/DBLE(num(1:r,2)-1))     ! SDev 
      tso(1:r)=crn(1:r,2)              ! Sort by index value
      CALL pair_sort(r,tso(1:r),tre(1:r)) 
      DO j=1,r
        p=tre(j) ; crn(j,4)=crn(p,3)   ! Sorted errors
        crn(p,5)=crn(p,3)/crn(p,2)     ! Not sorted, scaled errors
      ENDDO
      m=ad(nc)+yr(nc)-1
      CALL randnorm(m,ax)      ! Normally distributed numbers
      CALL quicksort(m,ax,ny)  ! Sort random
      CALL quicksort(m,dx,nx)  ! Sort tree indices 
      DO i=1,m ; tx(nx(i))=ax(ny(i)) ; ENDDO 
      crn(1:r,9:12)=0.D0 ; r=cyr(1)
      DO i=1,nc
        v=ly(i)-cfy(1)+1 ; u=MAX(1,v-yr(i)+1) 
        j=v-u+1 ; q=ad(i)+yr(i)-1 ; p=q-j+1
        WHERE (xok(p:q))
          crn(u:v,11)=crn(u:v,11)+tx(p:q)
          crn(u:v,12)=crn(u:v,12)+tx(p:q)**2
        END WHERE
      ENDDO    
      crn(1:r,9)=crn(1:r,11)/DBLE(num(1:r,2))      ! CRN
      crn(1:r,10)=SQRT(MAX(crn(1:r,12)-crn(1:r,9)* &  
        crn(1:r,11),0.001D0)/DBLE(num(1:r,2)-1))     ! SDev 
      WHERE (num(1:r,2).LT.4) crn(1:r,10)=0.D0
      tso(1:r)=crn(1:r,9)            ! Sort by index value
      CALL pair_sort(r,tso(1:r),tre(1:r)) 
      DO j=1,r ; crn(j,11)=crn(tre(j),10) ; ENDDO  ! Sorted normal
      OPEN(19,FILE="R2fig2.col",IOSTAT=ios,STATUS="REPLACE")
      WRITE(19,'("Column  1 = Year")')
      WRITE(19,'("Column  2 = Sample counts")')
      WRITE(19,'("Column  3 = One-curve RCS Chronology")')
      WRITE(19,'("Column  4 = S.Dev Sorted by Calendar Year")')
      WRITE(19,'("Column  5 = S.Dev Sorted on Index")')
      WRITE(19,'("Column  6 = S.Dev Scaled by Index")')
      WRITE(19,'("Column  7 = One-curve, Normal RCS Chronology")')
      WRITE(19,'("Column  8 = S.Dev Sorted by Calendar Year")')
      WRITE(19,'("Column  9 = S.Dev Sorted on Index")')
      DO i=1,cyr(1)
        WRITE(19,'(2I5,7F8.3)') cfy(1)-1+i,num(i,2),crn(i,2:5),crn(i,9:11)
      ENDDO
      CLOSE(19)
      RETURN
      END SUBROUTINE R2Fig2
!--------------------------------------------------------------
      SUBROUTINE RC_SM03()  ! Adjusted slope chronology 
      IMPLICIT NONE
      INTEGER              :: i,j,p,q,r,u,v
      REAL(8),DIMENSION(3) :: fac
      CALL det_default() ; idt=-2
      OPEN(74,FILE="zzz.prn",IOSTAT=ios,STATUS="REPLACE")
      fac=(/+0.001D0,0.D0,-0.001D0/) ; cf=1 
      DO j=1,3
        nc=0 ; CALL read_rft("../../raw/rcs/finnmrg.raw")
        CALL det_crnfy() ; r=cyr(1)
        crn(1:r,49+j)=(/(1.D0+DBLE(i)*fac(j),i=-r/2,r-r/2-1)/)
        DO i=1,nc
          r=yr(i) ; p=ad(i) ; q=p+r-1 
          u=fy(i)-cfy(1)+1 ; v=u+r-1
          x(p:q)=x(p:q)*crn(u:v,49+j)
        ENDDO
        IF (j.EQ.2) THEN
          cf=cf+5 ; idt=0 ; CALL detrend() ; cf=cf-5  ! cf=8
        ELSEIF (j.EQ.3) THEN
          cf=cf+4 ; idt=0 ; CALL detrend() ; cf=cf-4  ! cf=9
        ELSE
          idt=0 ; CALL detrend()         ! cf=1
        ENDIF
        idt=-2 ; cf=cf+1 ; sfo=1 ; CALL detrend()
        r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
        num(1:r,cf+10)=mcnt(1:r,mx)
        cf=cf+1 ; sfo=2 ; CALL detrend()
        r=sly(mx) ; crn(1:r,cf+10)=mval(1:r,mx)
        num(1:r,cf+10)=mcnt(1:r,mx)
        WRITE(74,'(4I6)') j,cfy(cf-2),cly(cf-2),cyr(cf-2)
        WRITE(74,'(4I6)') cf,cfy(cf-1),cly(cf-1),cyr(cf-1)
        WRITE(74,'(4I6)') cf,cfy(cf),cly(cf),cyr(cf)
      ENDDO
      r=cyr(1)
      crn(1:r,50)=crn(1:r,50)+SUM(crn(1:r,3)/crn(1:r,5))/DBLE(r) &
        -SUM(crn(1:r,50))/DBLE(r) 
      crn(1:r,52)=crn(1:r,52)+SUM(crn(1:r,7)/crn(1:r,5))/DBLE(r) &
        -SUM(crn(1:r,52))/DBLE(r) 
      DO i=1,12
        WRITE(74,'(4I6)') i,cfy(i),cly(i),cyr(i)
      ENDDO
      r=cyr(1)
      crn(1:r,36)=crn(1:r,6)/crn(1:r,4)
      crn(1:r,37)=crn(1:r,2)/crn(1:r,4)
      crn(1:r,38)=crn(1:r,7)/crn(1:r,5)
      crn(1:r,39)=crn(1:r,3)/crn(1:r,5)
      CLOSE(74)
      cnam(30)="Slope -0.001"
      cnam(31)="Slope +0.001"
      cnam(32)="Measured Values"
      cnam(33)="a) Finnish Trees - RCS Curves"
      cnam(34)="a) Finnish Trees - Ring-Width"
      RETURN
      END SUBROUTINE RC_SM03
!--------------------------------------------------------------
      SUBROUTINE RCS_RandSF(trg,yrs,sig) ! RCS SF using random trees 
      IMPLICIT NONE              
      INTEGER,INTENT(IN)     :: trg      ! Number of trees
      INTEGER,INTENT(IN)     :: yrs      ! Number of years
      INTEGER,INTENT(IN)     :: sig      ! Signal
      INTEGER                :: vnc      ! Stored common signal
      INTEGER,DIMENSION(mxs) :: vfy,vly,vyr,vad,vpth  ! Stored values
      REAL(8),DIMENSION(mxs) :: vpthr    ! Stored pith radius
      REAL(8),DIMENSION(mxd) :: vx       ! Stored signal free measures
      REAL(8),DIMENSION(mxs) :: wk       ! Stored pith radius
      INTEGER,DIMENSION(mxs) :: wki
      INTEGER                :: i,j,m,p,q
      REAL(8)                :: step,tot
      m=ad(nc)+yr(nc)-1      ; vx(1:m)=fx(1:m) ; vnc=nc ! Store signal free data
      vfy(1:nc)=fy(1:nc)     ; vly(1:nc)=ly(1:nc)
      vyr(1:nc)=yr(1:nc)     ; vpth(1:nc)=pth(1:nc)
      vpthr(1:nc)=pthr(1:nc) ; vad(1:nc+1)=ad(1:nc+1)
      step=DBLE(trg)/DBLE(yrs)         ! Step between each tree
      CALL RANDOM_NUMBER(wk(1:trg+1))  ! Select with replacement
      wki(1:trg+1)=NINT(wk(1:trg+1)*DBLE(vnc-1))+1
      tot=1.D0 ; nc=0
      DO j=1,yrs-60
        tot=tot+step ; m=INT(tot) ; tot=tot-DBLE(m)
        DO m=1,m       ! Random trees this year
          nc=nc+1 ; i=wki(nc) ; yr(nc)=vyr(i)
          fy(nc)=MIN(j,yrs-yr(nc)+1) ; ly(nc)=fy(nc)+yr(nc)-1 
          pth(nc)=fy(nc)-vfy(i)+vpth(i)
          p=ad(nc) ; q=p+yr(nc)-1 ; xok(p:q)=TR
          ad(nc+1)=q+1 ; pthr(nc)=vpthr(i)
          x(p:q)=vx(vad(i):vad(i)+yr(nc)-1)*crn(fy(nc):ly(nc),sig)
        ENDDO
      ENDDO 
      RETURN
      END SUBROUTINE RCS_RandSF
!--------------------------------------------------------------
      SUBROUTINE R2Fig8()  ! RCS Random chronologies 
      IMPLICIT NONE
      INTEGER,PARAMETER     :: trg=1200  ! Number of trees
      INTEGER,PARAMETER     :: yrs=2000  ! Number of years
      INTEGER,DIMENSION(12) :: seed
      INTEGER               :: i
      REAL(8)               :: ra,rb
      seed=101 ; CALL random_seed(put=seed)  ! Initialise random sequence
      CALL det_default() ; idt=-2 ; sfo=2
      src=1 ; srcno=4 ; idb=2  ! Normal
      nc=0 ; CALL read_rft("../../raw/rcs/yml-all.raw")
!     nc=0 ; CALL read_rft("../../raw/rcs/torn-all.raw")
      cf=2 ; CALL detrend() 
      DO i=1,yrs          ! Added signal
        crn(i,99)=DSIN(DBLE(i+750)*pi/500.D0)/2.D0+1.D0
      ENDDO
      crn(1:yrs,100)=crn(1:yrs,99)
      WHERE (crn(1:yrs,99).LT.1.D0) &
        crn(1:yrs,99)=1.D0/(2.D0-crn(1:yrs,99))
      gtr=1  ; cnam(20)="Ratios Sine Wave with Mean 1.0"
      CALL RCS_randSF(trg,yrs,99)
      src=1 ; srcno=1 ; cf=9  ; CALL detrend()
      CALL splinet(yrs,crn(1:yrs,9),500,crn(1:yrs,10))    
      src=2 ; srcno=2 ; cf=11 ; CALL detrend()
      CALL splinet(yrs,crn(1:yrs,11),500,crn(1:yrs,12))    
      src=2 ; srcno=3 ; cf=13 ; CALL detrend()
      CALL splinet(yrs,crn(1:yrs,13),500,crn(1:yrs,14))    
      gtr=2  ; cnam(21)="Ratios Sine Wave with Mean Single"
      src=1 ; srcno=1 ; cf=19  ; CALL detrend()
      CALL splinet(yrs,crn(1:yrs,19),500,crn(1:yrs,20))    
      src=2 ; srcno=2 ; cf=21 ; CALL detrend()
      CALL splinet(yrs,crn(1:yrs,21),500,crn(1:yrs,22))    
      src=2 ; srcno=3 ; cf=23 ; CALL detrend()
      CALL splinet(yrs,crn(1:yrs,23),500,crn(1:yrs,24))    
      ra=SUM(ABS(crn(1:yrs,100)-1.D0))
      rb=SUM(ABS(crn(1:yrs,10)))
      crn(1:yrs,101)=(crn(1:yrs,100)-1.D0)*rb/ra  ! Rescale signal
      OPEN(19,FILE="R2fig8.col",IOSTAT=ios,STATUS="REPLACE")
      WRITE(19,'("Column  1 = Calendar year")')
      WRITE(19,'("Column  2 = Sample count")')
      WRITE(19,'("Column  3 = Artificial sine wave")')
      WRITE(19,'("Column  3 = One-curve RCS")')
      WRITE(19,'("Column  3 = Two-curve RCS")')
      WRITE(19,'("Column  3 = Three-curve RCS")')
      WRITE(19,'("Column  3 = Means =1.0 One-curve RCS")')
      WRITE(19,'("Column  3 = Means =1.0 Two-curve RCS")')
      WRITE(19,'("Column  3 = Means =1.0 Three-curve RCS")')
      WRITE(19,'("Column  3 = Mean single One-curve RCS")')
      WRITE(19,'("Column  3 = Mean single Two-curve RCS")')
      WRITE(19,'("Column  3 = Mean single Three-curve RCS")')
      DO i=1,cyr(9)
        WRITE(19,'(2I5,10F8.3)') cfy(9)-1+i,num(i,1),crn(i,101), &
          crn(i,9),crn(i,11),crn(i,13),crn(i,10), &
          crn(i,12),crn(i,14),crn(i,20),crn(i,22),crn(i,24)
      ENDDO
      CLOSE(19)
      RETURN
      END SUBROUTINE R2Fig8
!--------------------------------------------------------------
      SUBROUTINE R2Fig8d()  ! Plots Best Fit
      IMPLICIT NONE                 
      REAL(8),DIMENSION(mxy) :: wka
      REAL(8)                :: ra,rb
      INTEGER                :: i=1,p,q,r
      ra=-0.8D0 ; rb=+0.8D0
      p=cfy(9) ; q=cly(9) ; r=cyr(9) 
      wka(1:r)=(/(DBLE(i),i=p,q)/)   
      grl=200 ; grr=2000 ; grt=140 ; grb=480
      CALL TICKS(2,'X')      ! No X ticks
      CALL LABELS('NONE','X')
      CALL NAME('','X')            ! Axis name
      CALL NAME('Index Value','Y') ! Axis name
      CALL plot_trees(r,num(1:r,9))  
      CALL tombox(p,q,ra,rb)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(red) ; CALL CURVE(wka(1:r),crn(1:r,101),r)
      CALL MESSAG("Added Sine wave",grl+1300,grt-45)
      CALL SETCLR(black) ; CALL CURVE(wka(1:r),crn(1:r,9),r)
      CALL MESSAG("a) One_curve RCS",grl+800,grt+30)
      CALL MESSAG(cnam(20),grl+200,grt-45)
      CALL ENDGRF()
      grt=490 ; grb=830
      CALL tombox(p,q,ra,rb)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(red) ; CALL CURVE(wka(1:r),crn(1:r,101),r)
      CALL SETCLR(black) ; CALL CURVE(wka(1:r),crn(1:r,11),r)
      CALL MESSAG("b) Two_curve RCS",grl+800,grt+30)
      CALL ENDGRF()
      grt=840 ; grb=1180
      CALL tombox(p,q,ra,rb)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(red) ; CALL CURVE(wka(1:r),crn(1:r,101),r)
      CALL SETCLR(black) ; CALL CURVE(wka(1:r),crn(1:r,13),r)
      CALL MESSAG("c) Three_curve RCS",grl+800,grt+30)
      CALL ENDGRF()

      grt=1190 ; grb=1530
      CALL tombox(p,q,ra,rb)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(red) ; CALL CURVE(wka(1:r),crn(1:r,101),r)
      CALL MESSAG("d) Means = 1.0",grl+30,grt+30)
      CALL SETCLR(grey) ; CALL CURVE(wka(1:r),crn(1:r,10),r)
      CALL MESSAG("One-curve",grl+650,grt+30)
      CALL SETCLR(blue) ; CALL CURVE(wka(1:r),crn(1:r,12),r)
      CALL MESSAG("Two-curve",grl+1050,grt+30)
      CALL SETCLR(green) ; CALL CURVE(wka(1:r),crn(1:r,14),r)
      CALL MESSAG("Three-curve",grl+1450,grt+30)
      CALL ENDGRF()

      grt=1540 ; grb=1880
      CALL LABELS('FLOAT','X')
      CALL NAME('Calendar Year','X')  ! Axis name
      CALL NAME('Ring Width','Y')     ! Axis name
      CALL tombox(p,q,ra,rb)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! GRIDLINES
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(red) ; CALL CURVE(wka(1:r),crn(1:r,101),r)
      CALL MESSAG("e) Mean Single",grl+30,grt+30)
      CALL SETCLR(grey) ; CALL CURVE(wka(1:r),crn(1:r,20),r)
      CALL MESSAG("One-curve",grl+650,grt+30)
      CALL SETCLR(blue) ; CALL CURVE(wka(1:r),crn(1:r,22),r)
      CALL MESSAG("Two-curve",grl+1050,grt+30)
      CALL SETCLR(green) ; CALL CURVE(wka(1:r),crn(1:r,24),r)
      CALL MESSAG("Three-curve",grl+1450,grt+30)
      CALL TICKS(10,'X')      ! No X ticks
      CALL ENDGRF()
      RETURN 
      END SUBROUTINE R2Fig8d
!-------------------------------------------------------------------
      SUBROUTINE R2Fig4d()  
      IMPLICIT NONE                 
      INTEGER,PARAMETER      :: leng=20    ! number of years
      INTEGER,PARAMETER      :: start=700  ! first year
      INTEGER,DIMENSION(8)   :: col
      REAL(8),DIMENSION(mxy) :: wka
      INTEGER                :: i=1,p,q,r,u,v,tno
      tno=nc
      col=(/blue,red,cyan,green,brown,black,yellow,grey/)
      wka(1:leng)=(/(DBLE(i),i=start,start+leng-1)/) 
      grl=1720 ; grr=2400 ; grt=140 ; grb=740
      CALL TICKS(5,'X')      ! No X ticks
      CALL NAME('Year','X')  ! Axis name
      CALL NAME('','Y')      ! Axis name
      CALL tombox(start,start+leng,-2.4D0,3.8D0)
      CALL SETCLR(grey) ; CALL GRID(1,1) ! Gridlines
      CALL LINWID(1) ; CALL HEIGHT(22)
      DO i=1,tno
        CALL SETCLR(col(MOD(i,8)+1))
        CALL CURVE(wka(2:leng),crn(2:leng,59+i),leng-1)
      ENDDO
      CALL SETCLR(black)
      CALL MESSAG("c) 50-year Normalised RCS Indices",grl-100,grt-40)  
      CALL ENDGRF() ; CALL LINWID(1) 
      grl=200 ; grr=880
      CALL NAME('Index Values','Y')   ! Axis name
      CALL LABELS('FLOAT','Y') 
      CALL tombox(start,start+leng,0.D0,3.D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! Gridlines
      CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(black)
      DO i=1,tno
        CALL SETCLR(col(MOD(i,8)+1))
        CALL CURVE(wka(2:leng),crn(2:leng,59+tno+i),leng-1)
      ENDDO
      CALL SETCLR(black)
      CALL MESSAG("a) 50-year Spline Indices",grl+15,grt-40)  
      CALL ENDGRF() ; CALL LINWID(1) 

      grl=900 ; grr=1580
      CALL NAME('','Y')   ! Axis name
      CALL LABELS('NONE','Y') 
      CALL tombox(start,start+leng,0.D0,3.D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! Gridlines
      CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(black)
      DO i=1,tno
        CALL SETCLR(col(MOD(i,8)+1))
        CALL CURVE(wka(2:leng),crn(2:leng,59+2*tno+i),leng-1)
      ENDDO
      CALL SETCLR(black)
      CALL MESSAG("b) RCS Indices",grl+15,grt-40)  
      CALL ENDGRF() ; CALL LINWID(1) 
      CALL LABELS('FLOAT','Y') 

      p=cfy(2) ; q=cly(2)-1 ; r=q-p+1
      u=cfy(1)+p-1 ; v=cfy(1)+q-1
      wka(1:r)=(/(DBLE(i),i=u,v)/) 
      grl=200 ; grr=2400 ; grt=840 ; grb=1140
      CALL LABELS('NONE','X')
      CALL NAME('','X')              ! Axis name
      CALL NAME('RBar Values','Y')   ! Axis name
      CALL tombox(u,v,0.1D0,0.72D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! Gridlines
      CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(black)
      CALL line_miss(r,wka(1:r),crn(p:q,8),okc(p:q,8))
      CALL MESSAG("d) Spline RBar",grl+700,grt+30)
      CALL SETCLR(red)
      CALL line_miss(r,wka(1:r),crn(p:q,11),okc(p:q,8))
      CALL MESSAG("RCS RBar",grl+1400,grt+30)
      CALL SETCLR(black)
      CALL MESSAG(cnam(60),grl+300,grt-45)
      CALL MESSAG(cnam(61),grl+1200,grt-45)
      CALL ENDGRF() 
      grt=1150 ; grb=1750
      CALL LABELS('FLOAT','X')
      CALL NAME('Calendar Year','X')               ! Axis name
      CALL NAME('EPS Values','Y')   ! Axis name
      CALL plot_trees(r,num(p:q,1))  
      CALL tombox(u,v,0.7D0,1.0D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! Gridlines
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(black)
      CALL line_miss(r,wka(1:r),crn(p:q,16),okc(p:q,8))
      CALL MESSAG("e) Spline EPS",grl+600,grb-45)
      CALL SETCLR(red) 
      CALL line_miss(r,wka(1:r),crn(p:q,17),okc(p:q,8))
      CALL MESSAG("RCS EPS",grl+1000,grb-45)
      CALL SETCLR(blue)
      CALL line_miss(r,wka(1:r),crn(p:q,19),okc(p:q,8))
      CALL MESSAG("RCS Adjusted EPS",grl+1800,grb-45)
      CALL TICKS(10,'X')     ! No X ticks
      CALL SETCLR(black) ; CALL ENDGRF() 
      RETURN 
      END SUBROUTINE R2Fig4d
!-------------------------------------------------------------------
      SUBROUTINE R2Fig4() ! Tornetrask TRW EPS 
      IMPLICIT NONE                 
      REAL(8),DIMENSION(mxd) :: zx
      INTEGER,DIMENSION(30)  :: tlist      ! tree list
      REAL(8)                :: mn,sd,sd1,sd2,cor,rb,eps
      INTEGER,PARAMETER      :: leng=50    ! number of years
      INTEGER,PARAMETER      :: start=700 ! first year
      INTEGER                :: i,j,p,q,tno,nr
      cf=5 ; CALL det_default() ; sfo=2 ; idt=-2 
      cnam(cf)="../../raw/rcs/torn-all.raw"
      wnam(cf)="Tornetrask TRW" 
      nc=0 ; CALL read_rft(cnam(cf))    ! Read raw data
      cnam(6)="EPS_TornTRW.prn"      
      OPEN(71,FILE=TRIM(cnam(6)),IOSTAT=ios,STATUS="REPLACE")
      src=2 ; srcno=2 ; idb=1    ! Two-curve Ratio RCS 
      cnam(7)="Two-Curve, Signal-Free, Ratio RCS"      
      cf=1 ; CALL EPS_prep1(50)    ! Process data
      cnam(49)="../../raw/rcs/yamalad.raw"
      WRITE(71,*)
      nc=0 ; CALL det_default() ; idt=-2 
      CALL read_rft(cnam(49)) 
      cf=51 ; idt=-2 ; CALL detrend()  ! RCS=1
      tno=0 ; p=start-cfy(51)+1 ; q=p+leng-1
      sd1=SUM(xcsd(p:q,mx))/DBLE(leng)
      i=ad(nc)+yr(nc)-1 ; zx(1:i)=dx(1:i)
      idt=52 ; cf=52 ; CALL detrend()  ! Spline=2
      sd2=SUM(xcsd(p:q,mx))/DBLE(leng)
      DO j=1,nc     ! Select trees to use
        IF (fy(j).LT.start.AND.ly(j).GT.start+leng) THEN
          tno=tno+1 ; tlist(tno)=j
        ENDIF
      ENDDO        
      DO i=1,tno
        j=tlist(i) ; p=ad(j)+start-fy(j) ; q=p+leng-1
        mn=SUM(dx(p:q))/DBLE(leng)
        sd=SQRT(SUM((dx(p:q)-mn)**2)/DBLE(leng-1))
        crn(1:leng,59+i)=(dx(p:q)-mn)/sd
        crn(1:leng,59+tno+i)=dx(p:q)
        crn(1:leng,59+2*tno+i)=zx(p:q)
      ENDDO
      rb=0.D0 ; nr=0    ! Rbar for normalised data
      DO i=1,tno-1 ; DO j=i+1,tno
        cor=SUM(crn(1:leng,59+i)*crn(1:leng,59+j))/ &
          (SQRT(SUM(crn(1:leng,59+i)**2)*SUM(crn(1:leng,59+j)**2)))
        rb=rb+cor ; nr=nr+1
      ENDDO ; ENDDO
      rb=rb/DBLE(nr)    ! EPS=(N*RB)/(1+(N-1)*RB
      eps=DBLE(tno)*rb/(1.D0+DBLE(tno-1)*rb) ! EPS & Divide by root N#
      mn=(1.D0-rb)*0.9D0/(rb*(1.D0-0.9D0))
      WRITE(71,'("start, tno, sd1, sd2, (sd1/sd2)**2")') 
      WRITE(71,'("rb, eps, mn, mn*(sd1/sd2)**2")') 
      WRITE(71,'(2I6,5F8.3,2F8.1)')  &
         start,tno,sd1,sd2,(sd1/sd2)**2,rb,eps,mn,mn*(sd1/sd2)**2
      CLOSE(71) ; nc=tno
      OPEN(19,FILE="R2fig4.col",IOSTAT=ios,STATUS="REPLACE")
      WRITE(19,'("Panel a)")')
      WRITE(19,'("Column  1 = Year")')
      WRITE(19,'("Columns 2+ are each tree")')
      DO i=2,leng
        WRITE(19,'(I5,12F8.3)') start-1+i,crn(i,60:59+tno)
      ENDDO
      WRITE(19,'("Panel b)")')
      WRITE(19,'("Column  1 = Year")')
      WRITE(19,'("Columns 2+ are each tree")')
      DO i=2,leng
        WRITE(19,'(I5,12F8.3)') start-1+i,crn(i,60+tno:59+2*tno)
      ENDDO
      WRITE(19,'("Panel b)")')
      WRITE(19,'("Column  1 = Year")')
      WRITE(19,'("Columns 2+ are each tree")')
      DO i=2,leng
        WRITE(19,'(I5,12F8.3)') start-1+i,crn(i,60+2*tno:59+3*tno)
      ENDDO
      WRITE(19,'("Column  1 = Year")')
      WRITE(19,'("Column  2 = Sample counts")')
      WRITE(19,'("Column  3 = Spline RBar")')
      WRITE(19,'("Column  4 = RCS RBar")')
      WRITE(19,'("Column  5 = Spline EPS")')
      WRITE(19,'("Column  6 = RCS EPS")')
      WRITE(19,'("Column  7 = RCS Adjusted EPS")')
      DO i=cfy(2),cly(2)
        WRITE(19,'(2I5,7F8.3)') cfy(1)-1+i,num(i,1),crn(i,8), &
           crn(i,11),crn(i,16),crn(i,17),crn(i,18)
      ENDDO
      CLOSE(19)
      RETURN 
      END SUBROUTINE R2Fig4
!------------------------------------------------------------------------
      SUBROUTINE R2Fig1d()   ! Plots Pith/not paper figure 1
      IMPLICIT NONE                 
      REAL(8),DIMENSION(mxy) :: wka
      REAL(8)                :: ra,rb
      INTEGER                :: i=1,j=1,p,q,r
      grl=200 ; grr=2000 ; grt=140 ; grb=400
      CALL NAME('Ring Width','Y') ! Axis name
      CALL NAME('','X')           ! Axis name
      CALL LABELS('NONE','X')     ! Axis labels
      DO j=1,3
        IF     (j.EQ.1) THEN
          ra=0.5D0 ; rb=0.7D0 ; r=450
          CALL TICKS(5,'X')          ! X ticks 
        ELSEIF (j.EQ.2) THEN
          ra=0.2D0 ; rb=1.5D0 ; r=450
          CALL LABELS('FLOAT','X')   ! Axis labels
        ELSEIF (j.EQ.3) THEN
          ra=0.0D0 ; rb=0.7D0 ; r=1045
          CALL NAME('Ring Age','X')  ! Axis name
          CALL TICKS(2,'Y')          ! Y ticks 
          CALL TICKS(10,'X')         ! X ticks 
        ENDIF
        wka(1:r)=(/(DBLE(i),i=1,r)/) 
        CALL plot_trees(r,num(1:r,j+9))  
        CALL tombox(1,r,ra,rb)
        CALL SETCLR(grey) ; CALL GRID(1,1)  ! GRIDLINES
        CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(blue)
        p=cfy(j+3) ; q=r ; r=q-p+1
        CALL CURVE(wka(p:q),crn(p:q,j+3),r)
        IF (j.EQ.1) CALL MESSAG("Without Pith",grl+700,grt-40)
        CALL SETCLR(red) ; CALL CURVE(wka(p:q),crn(p:q,j+9),r)
        IF (j.EQ.1) CALL MESSAG("With Pith",grl+1200,grt-40)
        CALL SETCLR(black) ; CALL MESSAG(cnam(j+9),grl+750,grt+30)
        CALL ENDGRF() 
        IF (j.EQ.2) THEN
          grt=grt+320 ; grb=grb+320
        ELSE
          grt=grt+270 ; grb=grb+270
        ENDIF
      ENDDO
      grt=grt+120 ; grb=grb+120
      CALL NAME('','X')        ! Axis name
      CALL LABELS('NONE','X')  ! Axis labels
      DO j=1,3
        r=cyr(j) ; wka(1:r)=(/(DBLE(i),i=cfy(j),cly(j))/)  
        CALL plot_trees(r,num(1:r,j))  
        CALL NAME('Index','Y')   ! Axis name
        IF     (j.EQ.1) THEN
          CALL TICKS(4,'X')         ! X ticks 
          ra=0.8D0 ; rb=1.2D0
        ELSEIF (j.EQ.2) THEN
          ra=0.4D0 ; rb=1.8D0
          CALL LABELS('FLOAT','X')  ! Axis labels
          CALL TICKS(5,'Y')         ! Y ticks 
        ELSEIF (j.EQ.3) THEN
          ra=0.4D0 ; rb=1.99D0
          CALL NAME('Calendar Year','X')  ! Axis name
        ENDIF
        CALL tombox(cfy(j),cly(j),ra,rb)
        CALL SETCLR(grey) ; CALL GRID(1,1)  ! GRIDLINES
        CALL LINWID(1) ; CALL HEIGHT(22)
        CALL SETCLR(blue)  ; CALL CURVE(wka(1:r),crn(1:r,j),r)
        CALL SETCLR(red)   ; CALL CURVE(wka(1:r),crn(1:r,6+j),r)
        CALL SETCLR(black) ; CALL MESSAG(cnam(j+12),grl+800,grt+30)
        CALL ENDGRF()  
        IF (j.EQ.2) THEN
          grt=grt+320 ; grb=grb+320
        ELSE
          grt=grt+270 ; grb=grb+270
        ENDIF
      ENDDO
      CALL TICKS(10,'X')         ! X ticks 
      RETURN
      END SUBROUTINE R2Fig1d
!-------------------------------------------------------------------
      SUBROUTINE R2Fig1() 
      IMPLICIT NONE                 
      INTEGER            :: i,j,m,n,r
      cnam(40)="../../raw/rcs/S88G1112A.mxd"
      cnam(10)="a) Tornetrask MXD RCS Curves"
      cnam(13)="d) Tornetrask MXD Chronology"
      cnam(41)="../../raw/rcs/S88G0812.raw"
      cnam(11)="b) Tornetrask TRW RCS Curves"
      cnam(14)="e) Tornetrask TRW Chronology"
      cnam(42)="../../raw/rcs/chin005c.rwl"
      cnam(12)="c) Chin005 TRW RCS Curves"
      cnam(15)="f) Chin005 TRW Chronology"
      CALL det_default() ; idt=-2 ; CDsp=30
      DO j=1,3
        nc=0 ; cf=j ; CALL read_rft(cnam(39+j))
        poo=2 ; CALL detrend()  ! Pith offset OFF
        r=cyr(cf) ; crn(1:r,cf)=xcsm(1:r,mx)
        cf=cf+3 ; r=sly(mx) ; crn(1:r,cf)=msmo(1:r,mx)
        num(1:r,cf)=mcnt(1:r,mx)
        cfy(cf)=sfy(mx) ; cly(cf)=r ; cyr(cf)=r
        poo=1 ; cf=cf+3 ; CALL detrend()  ! Pith offset ON
        r=cyr(cf) ; crn(1:r,cf)=xcsm(1:r,mx)
        cf=cf+3 ; r=sly(mx) ; crn(1:r,cf)=msmo(1:r,mx)
        num(1:r,cf)=mcnt(1:r,mx)
        cfy(cf)=sfy(mx) ; cly(cf)=r ; cyr(cf)=r
      ENDDO   
      OPEN(19,FILE="R2fig1.col",IOSTAT=ios,STATUS="REPLACE")
      WRITE(19,'("Column  1 = Ring Age")')
      WRITE(19,'("Column  2 = Chin005 TRW sample count")')
      WRITE(19,'("Column  3 = Chin005 TRW No Pith")')
      WRITE(19,'("Column  4 = Chin005 TRW With Pith")')
      WRITE(19,'("Column  5 = Torn MXD sample count")')
      WRITE(19,'("Column  6 = Torn MXD No Pith")')
      WRITE(19,'("Column  7 = Torn MXD With Pith")')
      WRITE(19,'("Column  8 = Torn TRW sample count")')
      WRITE(19,'("Column  9 = Torn TRW No Pith")')
      WRITE(19,'("Column 10 = Torn TRW With Pith")')
      DO i=1,cyr(4)
        WRITE(19,'(I4,3(I8,2F8.3))') i,num(i,6),crn(i,6),crn(i,12), &
          num(i,4),crn(i,4),crn(i,10),num(i,5),crn(i,5),crn(i,11)
      ENDDO
      DO i=cyr(4)+1,cyr(6)
        WRITE(19,'(I4,I8,2F8.3)') i,num(i,6),crn(i,6),crn(i,12)
      ENDDO
      WRITE(19,*)
      WRITE(19,'("Column  1 = Year")')
      WRITE(19,'("Column  2 = Chin005 TRW sample count")')
      WRITE(19,'("Column  3 = Chin005 TRW No Pith")')
      WRITE(19,'("Column  4 = Chin005 TRW With Pith")')
      WRITE(19,'("Column  5 = Torn MXD sample count")')
      WRITE(19,'("Column  6 = Torn MXD No Pith")')
      WRITE(19,'("Column  7 = Torn MXD With Pith")')
      WRITE(19,'("Column  8 = Torn TRW sample count")')
      WRITE(19,'("Column  9 = Torn TRW No Pith")')
      WRITE(19,'("Column 10 = Torn TRW With Pith")')
      DO i=cfy(1),cfy(3)-1
        m=i-cfy(3)+1 ; n=i-cfy(1)+1
        WRITE(19,'(I4,24X,2(I8,2F8.3))') i, &
          num(n,7),crn(n,1),crn(n,7),num(n,8),crn(n,2),crn(n,8)
      ENDDO
      DO i=cfy(3),cly(3)
        m=i-cfy(3)+1 ; n=i-cfy(1)+1
        WRITE(19,'(I4,3(I8,2F8.3))') i,num(m,3),crn(m,3),crn(m,9), &
          num(n,1),crn(n,1),crn(n,7),num(n,2),crn(n,2),crn(n,8)
      ENDDO
      DO i=cly(3)+1,cly(1)
        m=i-cfy(3)+1 ; n=i-cfy(1)+1
        WRITE(19,'(I4,24X,2(I8,2F8.3))') i, &
          num(n,1),crn(n,1),crn(n,7),num(n,2),crn(n,2),crn(n,8)
      ENDDO
      CLOSE(19)
      RETURN 
      END SUBROUTINE R2Fig1
!------------------------------------------------------------------------
      SUBROUTINE smerr1() ! Slope, Mean and error (normal distribution)
      IMPLICIT NONE                 
      REAL(8),DIMENSION(mxt) :: wka
      REAL(8)                :: rr,mn,sl,yi,xi 
      INTEGER                :: i=1,p,q,r,u,v
!     sd=SQRT(SUM((dx(1:m)-mn)**2)/(rr-1.D0))
      r=cyr(cf) ; crn(1:r,cf+1:cf+10)=0.D0 
      DO i=1,nc     ! For each tree
        p=ad(i) ; r=yr(i) ; q=p+r-1      ! Ring address
        u=fy(i)-cfy(cf)+1 ; v=u+r-1      ! Chronology address
        rr=DBLE(COUNT(xok(p:q)))
        mn=SUM(dx(p:q),MASK=xok(p:q))/rr ! Mean  indices
        ax(p:q)=dx(p:q)-mn               ! Remove mean
        CALL TRENDmiss(r,ax(p:q),wka(1:r),xok(p:q),sl,yi,xi)  ! Best fit line
        WHERE (xok(p:q))
          crn(u:v,cf+4)=crn(u:v,cf+4)+dx(p:q)      ! Sum of values    
          crn(u:v,cf+5)=crn(u:v,cf+5)+mn           ! Sum of means
          crn(u:v,cf+6)=crn(u:v,cf+6)+wka(1:r)     ! Sum of slopes 
          crn(u:v,cf+7)=crn(u:v,cf+7)+dx(p:q)**2   ! Sums of squares
          crn(u:v,cf+8)=crn(u:v,cf+8)+mn**2
          crn(u:v,cf+9)=crn(u:v,cf+9)+wka(1:r)**2
        END WHERE 
      ENDDO  
      r=cyr(cf)  
      crn(1:r,cf+10)=DBLE(num(1:r,cf))    ! Floating counts
      WHERE (num(1:r,cf).GE.1)
        crn(1:r,cf+1)=crn(1:r,cf+4)/crn(1:r,cf+10)  ! Mean
        crn(1:r,cf+2)=crn(1:r,cf+5)/crn(1:r,cf+10)  ! Mean
        crn(1:r,cf+3)=crn(1:r,cf+6)/crn(1:r,cf+10)  ! Mean
      END WHERE
      crn(1:r,cf+10)=crn(1:r,cf+10)-1.D0   ! Counts -1
      WHERE (num(1:r,cf).GT.3) 
        crn(1:r,cf+4)=SQRT(MAX(crn(1:r,cf+7)-crn(1:r,cf+1)* &  
          crn(1:r,cf+4),0.001D0)/crn(1:r,cf+10))        ! SDev 
        crn(1:r,cf+5)=SQRT(MAX(crn(1:r,cf+8)-crn(1:r,cf+2)* &  
          crn(1:r,cf+5),0.001D0)/crn(1:r,cf+10))        ! SDev 
        crn(1:r,cf+6)=SQRT(MAX(crn(1:r,cf+9)-crn(1:r,cf+3)* &  
          crn(1:r,cf+6),0.001D0)/crn(1:r,cf+10))        ! SDev 
      ELSEWHERE 
        crn(1:r,cf+4)=0.D0 ; crn(1:r,cf+5)=0.D0 ; crn(1:r,cf+6)=0.D0 
      END WHERE    
      DO i=cf+1,cf+6    ! 100-year smoothing
        CALL spline_miss(r,crn(1:r,i),50,crn(1:r,i),num(1:r,cf).GE.1) 
      ENDDO
      RETURN 
      END SUBROUTINE smerr1
!------------------------------------------------------------------------
      SUBROUTINE smerrd(ref1)  
      IMPLICIT NONE                 
      INTEGER,INTENT(IN)     :: ref1
      REAL(8),DIMENSION(mxy) :: wka 
      INTEGER                :: i=1,j,k,m,p,q,r
      REAL(8)                :: ra,rc
      INTEGER,DIMENSION(4)   :: col
      col=(/blue,red,cyan,black/)
      grl=200 ; grr=2000 ; grt=150 ; grb=530
      CALL TICKS(5,'X')        ! X ticks 
      CALL NAME('','X')
      CALL LABELS('NONE','X')
      q=cyr(1) ; wka(1:q)=(/(DBLE(i),i=cfy(1),cly(1))/) 
      IF (ref1.LT.3) THEN  ! Ignore early low counts
        p=1-cfy(1)+1 ; j=ref1*3-2
      ELSE 
        p=-1500-cfy(1)+1 ; j=ref1*3-8
      ENDIF
      r=q-p+1
      DO i=j,j+2
        IF (i.EQ.1.OR.i.EQ.4) &
          CALL plot_trees(r,num(p:q,1))  
        IF (i.EQ.3.OR.i.EQ.6) THEN
          CALL NAME('Calendar Year','X')
          CALL LABELS('FLOAT','X')
        ENDIF
        CALL NAME('Index','Y')   ! Axis name

        IF (ref1.EQ.1) THEN   
          IF (i.EQ.j) THEN   
            CALL tombox(1,cly(1),-0.95D0,1.05D0)
          ELSE
            CALL tombox(1,cly(1),-0.25D0,0.75D0)
          ENDIF
        ELSEIF (ref1.EQ.2) THEN
          IF (i.EQ.4) THEN   
            CALL tombox(1,cly(1),0.58D0,+1.0D0)
          ELSE
            CALL tombox(1,cly(1),0.12D0,0.65D0)
          ENDIF
        ELSEIF (ref1.EQ.3) THEN
          IF (i.EQ.1) THEN   
            CALL tombox(-1500,cly(1),-0.8D0,1.0D0)
          ELSE
            CALL tombox(-1500,cly(1),-0.38D0,0.3D0)
          ENDIF
        ELSEIF (ref1.EQ.4) THEN
          IF (i.EQ.4) THEN   
            CALL tombox(-1500,cly(1),0.5D0,+1.1D0)
          ELSE 
            CALL tombox(-1500,cly(1),0.18D0,0.69D0)
          ENDIF
        ENDIF
        CALL SETCLR(grey) ; CALL GRID(1,1)       ! Gridlines
        CALL LINWID(1) ; CALL HEIGHT(22)
        DO k=1,4
          CALL SETCLR(col(k))
          IF (i.EQ.1.OR.i.EQ.4) &
            CALL MESSAG(cnam(30+k),grl+150+k*350,grt-45)
          CALL CURVE(wka(p:q),crn(p:q,i+k*10-9),r)
        ENDDO 
        IF (i.EQ.1.OR.i.EQ.4) CALL MESSAG(cnam(30),grl+10,grt-45)
        CALL MESSAG(cnam(41+i-j),grl+700,grt+30)
        CALL ENDGRF() 
        grt=grt+390 ; grb=grb+390
        CALL TICKS(2,'Y')      ! Y ticks 
      ENDDO 
      CALL TICKS(10,'X')       ! Y ticks 
      CALL TICKS(5,'Y')        ! Y ticks 
      IF (ref1.EQ.2.OR.ref1.EQ.4) RETURN ! no frequency on S.Dev.
      grt=grt+130 ; grb=grb+200
      m=cyr(49)   ! Number of lags
      ra=LOG10(MINVAL(crn(2:m,52))) ; rc=DBLE(INT(ra))
      CALL NAME('Period (years)','X')
      CALL NAME('Power','Y')   ! Axis name
      CALL NEGLOG(0.00001D0)
      CALL LABELS('FLOAT','X')
      CALL AXSPOS(grl,grb)
      CALL AXSLEN(grr-grl,grb-grt)
      CALL HEIGHT(18) 
      CALL SETCLR(black)
      CALL AXSSCL('LOG','XY')   ! logscale XY
      CALL LABDIG(3,'Y')
      CALL GRAF(LOG10(5.D0),LOG10(DBLE(m+10)), &
        1.D0,1.D0,ra,LOG10(3.D0),rc,1.0D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)    ! Gridlines
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(black)
      m=(m*39)/100
      DO i=1,4
        j=44+i*6 ; CALL SETCLR(col(i))
        CALL curve(1.D0/crn(2:m,j+1),crn(2:m,j+2),m-1)
      ENDDO 
      CALL SETCLR(black)
      CALL MESSAG("d) Power spectrum plot",grl+100,grt+30)
      CALL ENDGRF() 
      CALL AXSSCL('LIN','XY')  ! linear scale XY
      RETURN 
      END SUBROUTINE smerrd
!------------------------------------------------------------------------
      SUBROUTINE smerr(ref1)  ! Create chronologies + St Err
      IMPLICIT NONE                 
      CHARACTER(10),DIMENSION(4),PARAMETER :: fnam= &
        (/"R2fig6.col","R2sm2.col ","R2fig7.col","R2sm3.col "/)
      INTEGER,INTENT(IN) :: ref1
      INTEGER            :: i,j,m,p,q,r,lag
      IF (ref1.LT.3) THEN
        cnam(1)="../../raw/rcs/yml-all.raw"  ! All trees
      ELSE
        cnam(1)="../../raw/rcs/Allm.raw"  ! All trees
      ENDIF
      nc=0 ; CALL read_rft(cnam(1))
      CALL det_default() ; idt=-2 ; idb=2 ! Normal distribution
      cnam(30)="Normal SF RCS"
      cnam(31:34)= &
        (/"One-curve  ","Two-curve  ","Three-curve","Four-curve "/)
      src=1 ; srcno=1 ; cf=1
      DO i=1,4
        CALL detrend() ; CALL smerr1()      
        IF (i.EQ.1) THEN  ! First time
          src=2 
          PD: DO p=1,cyr(1)   ! Find > 5 trees section
            IF (num(p,1).GE.5) EXIT PD
          ENDDO PD
          q=cyr(1) ; r=q-p+1 ; lag=r/4
          cyr(49)=lag    ; cfy(50)=p-1+cfy(1)
          cly(50)=cly(1) ; cyr(50)=r
        ENDIF
        IF (ref1.EQ.1.OR.ref1.EQ.3) THEN ! no frequency on S.Dev.
          m=44+i*6 ; crn(1:r,m)=crn(p:q,cf)
          CALL frequency(r,lag,crn(1:r,m),crn(1:lag+1,m+1:m+5))
        ENDIF
        srcno=i+1 ; cf=cf+10
      ENDDO
      IF (ref1.EQ.1.OR.ref1.EQ.3) THEN 
        cnam(41)="a) 50-year Low Pass"
        cnam(42)="b) Means Only"
        cnam(43)="c) Slopes Only"
      ELSE
        cnam(41)="a) 50-year Low Pass Std. dev"
        cnam(42)="b) Means Only  Std. dev"
        cnam(43)="c) Slopes Only Std. dev"
      ENDIF
      r=cyr(1)
      OPEN(19,FILE=fnam(ref1),IOSTAT=ios,STATUS="REPLACE")
      WRITE(19,'("Column  1 = Calendar year")')
      WRITE(19,'("Column  2 = Sample count")')
      WRITE(19,'("Column  3 = ",A)') TRIM(cnam(41)(4:28))//" "//TRIM(cnam(31))
      WRITE(19,'("Column  4 = ",A)') TRIM(cnam(41)(4:28))//" "//TRIM(cnam(32))
      WRITE(19,'("Column  5 = ",A)') TRIM(cnam(41)(4:28))//" "//TRIM(cnam(33))
      WRITE(19,'("Column  6 = ",A)') TRIM(cnam(41)(4:28))//" "//TRIM(cnam(34))
      WRITE(19,'("Column  7 = ",A)') TRIM(cnam(42)(4:28))//" "//TRIM(cnam(31))
      WRITE(19,'("Column  8 = ",A)') TRIM(cnam(42)(4:28))//" "//TRIM(cnam(32))
      WRITE(19,'("Column  9 = ",A)') TRIM(cnam(42)(4:28))//" "//TRIM(cnam(33))
      WRITE(19,'("Column 10 = ",A)') TRIM(cnam(42)(4:28))//" "//TRIM(cnam(34))
      WRITE(19,'("Column 11 = ",A)') TRIM(cnam(43)(4:28))//" "//TRIM(cnam(31))
      WRITE(19,'("Column 12 = ",A)') TRIM(cnam(43)(4:28))//" "//TRIM(cnam(32))
      WRITE(19,'("Column 13 = ",A)') TRIM(cnam(43)(4:28))//" "//TRIM(cnam(33))
      WRITE(19,'("Column 14 = ",A)') TRIM(cnam(43)(4:28))//" "//TRIM(cnam(34))
      j=0 ; IF (ref1.EQ.2.OR.ref1.EQ.4) j=3 
      DO i=1,cyr(1)
        WRITE(19,'(2I5,12F8.3)') cfy(1)-1+i,num(i,1), &
          crn(i,2+j),crn(i,12+j),crn(i,22+j),crn(i,32+j), &
          crn(i,3+j),crn(i,13+j),crn(i,23+j),crn(i,33+j), &
          crn(i,4+j),crn(i,14+j),crn(i,24+j),crn(i,34+j)
      ENDDO
      IF (ref1.EQ.1.OR.ref1.EQ.3) THEN
        WRITE(19,'("Column  1 = Number")')
        WRITE(19,'("Column  2 = Period ",A)') 
        WRITE(19,'("Column  3 = Power  ",A)') TRIM(cnam(31))
        WRITE(19,'("Column  5 = Power  ",A)') TRIM(cnam(32))
        WRITE(19,'("Column  7 = Power  ",A)') TRIM(cnam(33))
        WRITE(19,'("Column  9 = Power  ",A)') TRIM(cnam(34))
        DO i=2,(cyr(49)*39)/100
          WRITE(19,'(I4,5F10.3)') i,1.D0/crn(i,51),crn(i,52),  &
            crn(i,58),crn(i,64),crn(i,70)
        ENDDO
      ENDIF
      CLOSE(19)
      RETURN 
      END SUBROUTINE smerr
!------------------------------------------------------------------------
      SUBROUTINE R2SM1d()  
      IMPLICIT NONE                 
      CHARACTER(26),DIMENSION(3:5) :: lab
      REAL(8),DIMENSION(mxy)       :: wka
      INTEGER :: i=1,j,r,ra
      lab=(/"b) S.Dev. by Calendar Year", &
            "c) S.Dev. Sorted on Index ", &
            "d) S.Dev. Scaled by Index "/)
      r=cyr(1) ; wka(1:r)=(/(DBLE(i),i=cfy(1),cly(1))/) 
      grl=200 ; grr=2400 ; ra=140 ; grt=ra ; grb=ra+240
      CALL plot_trees(r,num(1:r,2))  
      CALL LABELS('NONE','X')
      CALL TICKS(2,'Y') 
      CALL NAME('','X')        ! Axis name
      CALL NAME('Index','Y')   ! Axis name
      CALL tombox(cfy(1),cly(1),0.0D0,1.99D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)       ! Gridlines
      CALL LINWID(1) ; CALL HEIGHT(22) ; CALL SETCLR(black)
      CALL CURVE(wka(1:r),crn(1:r,2),r)
      CALL MESSAG(cnam(21),grl+900,grt-45)
      CALL MESSAG(wnam(21),grl+100,grt-45)
      CALL ENDGRF() ; CALL LINWID(1) 
      ra=ra+10
      DO j=3,4
        ra=ra+240 ; grt=ra ; grb=ra+230
        IF (j.EQ.4) THEN
          CALL LABELS('FLOAT','X')
          CALL NAME('Calendar Year','X')               ! Axis name
        ENDIF
        CALL tombox(cfy(1),cly(1),0.D0,1.5D0)
        CALL LINWID(4) ; CALL SETCLR(silver)
        DO i=1,r
          CALL RLINE(wka(i),0.D0,wka(i),crn(i,j))
        ENDDO
        CALL LINWID(1) ; CALL SETCLR(grey) ; CALL GRID(1,1)       ! Gridlines
        CALL HEIGHT(22) ; CALL SETCLR(black)
        CALL MESSAG(lab(j),grl+900,grt+30)
        CALL ENDGRF() 
      ENDDO
      CALL TICKS(5,'Y') 
      RETURN 
      END SUBROUTINE R2SM1d
!-------------------------------------------------------------------
      SUBROUTINE R2SM1()   ! yaqmal SDev power transform
      IMPLICIT NONE
      INTEGER  :: i,j,p,q,r,u,v
      cnam(1)="../../raw/rcs/yml-all.raw"
      wnam(21)="Yamal TRW"
      nc=0 ; CALL read_rft(cnam(1))
      CALL det_default() ; idt=-2 ; cf=1 
      src=1 ; srcno=1 ; itn=2 ; ind=2 
      cnam(21)="a) One-curve RCS with Power Transform"
      CALL detrend() ; r=cyr(1)
      crn(1:r,2:8)=0.D0 ; num(1:r,2)=0
      DO i=1,nc              ! Simple RCS
        p=ad(i) ; j=yr(i) ; q=p+j-1 ; u=fy(i)-cfy(1)+1 ; v=u+j-1
        WHERE (xok(p:q))
          num(u:v,2)=num(u:v,2)+1
          crn(u:v,3)=crn(u:v,3)+dx(p:q)
          crn(u:v,4)=crn(u:v,4)+dx(p:q)**2
        END WHERE
      ENDDO                  ! Ignore counts < 4
      JD: DO j=1,r ; IF (num(j,2).GE.4) EXIT JD ; ENDDO JD
      crn(1:r-j+1,1:4)=crn(j:r,1:4) ; num(1:r-j+1,2)=num(j:r,2)  
      cfy(1)=cfy(1)+j-1 ; r=r-j+1 ; cyr(1)=r
      crn(1:r,2)=crn(1:r,3)/DBLE(num(1:r,2))      ! CRN
      crn(1:r,3)=SQRT(MAX(crn(1:r,4)-crn(1:r,2)* &  
        crn(1:r,3),0.001D0)/DBLE(num(1:r,2)-1))   ! SDev 
      tso(1:r)=crn(1:r,2)              ! Sort by index value
      CALL pair_sort(r,tso(1:r),tre(1:r)) 
      DO j=1,r
        p=tre(j) ; crn(j,4)=crn(p,3)   ! Sorted errors
        crn(p,5)=crn(p,3)/crn(p,2)     ! Not sorted, scaled errors
      ENDDO
      tso(1:r)=DBLE(num(1:r,2))        ! Sort by sample count
      CALL pair_sort(r,tso(1:r),tre(1:r)) 
      DO j=1,r
        p=tre(j) ; crn(j,6)=crn(p,5)   ! Sorted errors
        crn(p,7)=crn(p,5)*SQRT(DBLE(num(p,2)))/SQRT(30.D0)
      ENDDO    ! Not sorted, scaled errors
      idb=1   ! Normal index distribution off
      OPEN(19,FILE="R2sm1.col",IOSTAT=ios,STATUS="REPLACE")
      WRITE(19,'("Column  1 = Year")')
      WRITE(19,'("Column  2 = Sample counts")')
      WRITE(19,'("Column  3 = One-curve RCS with Power Transform")')
      WRITE(19,'("Column  4 = S.Dev Sorted by Calendar Year")')
      WRITE(19,'("Column  5 = S.Dev Sorted on Index")')
      WRITE(19,'("Column  6 = S.Dev Scaled by Index")')
      DO i=1,cyr(1)
        WRITE(19,'(2I5,3F8.3)') cfy(1)-1+i,num(i,2),crn(i,2:4)
      ENDDO
      CLOSE(19)
      RETURN
      END SUBROUTINE R2SM1
!--------------------------------------------------------------
      SUBROUTINE R2Fig5d()  ! Briffa 2013 - SM Figure YT14
      IMPLICIT NONE                 
      REAL(8),DIMENSION(mxy) :: wka 
      INTEGER                :: i=1,p,q,r,w
      INTEGER,DIMENSION(4) :: col
      col=(/blue,red,cyan,black/)
      w=370 ; wka(1:w)=(/(DBLE(i),i=1,w)/) 
      grl=200 ; grr=2000 ; grt=130 ; grb=450
      CALL TICKS(5,'X')    ! No X ticks
      CALL NAME('Ring Age','X')
      CALL NAME('Ring Width','Y')    ! Axis name
      CALL plot_trees(w,mcnt(1:w,mx))  
      CALL tombox(1,w,0.D0,1.7D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)  ! Gridlines
      CALL LINWID(1) ; CALL HEIGHT(22)
      DO i=2,5
        CALL SETCLR(col(i-1))
        p=cfy(i) ; q=MIN(cly(i),w) ; r=q-p+1
        CALL thickthin(r,wka(p:q),crn(p:q,i),num(p:q,i),7)
        CALL MESSAG(wnam(i),grl+i*400-750,grt-45)
      ENDDO 
      CALL SETCLR(black)
      CALL MESSAG("a) Three-curve RCS, Sig-free",grl+900,grt+30)
      CALL ENDGRF() 

      r=cyr(1) ; wka(1:r)=(/(DBLE(i),i=cfy(1),cly(1))/) 
      grt=580 ; grb=900
      p=1-cfy(1)+1 ; q=cyr(1) ; r=q-p+1 
      CALL LABELS('NONE','X')
      CALL NAME('','X')
      CALL NAME('Index Value','Y')    ! Axis name
      CALL tombox(1,cly(1),-1.2D0,1.2D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)  ! Gridlines
      CALL LINWID(1) ; CALL HEIGHT(22)
      DO i=6,8
        CALL SETCLR(col(i-5))
        CALL thickthin(r,wka(p:q),crn(p:q,i),num(p:q,i),7)
      ENDDO 
      CALL SETCLR(black)
      CALL MESSAG("b) Chronology means 1.0",grl+900,grt+30)
      CALL ENDGRF() 

      grt=910 ; grb=1230
      p=1-cfy(1)+1 ; q=cyr(1) ; r=q-p+1 
      CALL NAME('Index Value','Y')    ! Axis name
      CALL tombox(1,cly(1),-1.2D0,1.2D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)  ! Gridlines
      CALL LINWID(1) ; CALL HEIGHT(22)
      DO i=6,8
        CALL SETCLR(col(i-5))
        CALL thickthin(r,wka(p:q),crn(p:q,i+10),num(p:q,i),7)
      ENDDO 
      CALL SETCLR(black)
      CALL MESSAG("c) Chronology Means Reset",grl+900,grt+30)
      CALL ENDGRF()

      r=cyr(1) ; wka(1:r)=(/(DBLE(i),i=cfy(1),cly(1))/) 
      grt=1240 ; grb=1560
      p=1-cfy(1)+1 ; q=cyr(1) ; r=q-p+1 
      CALL LABELS('NONE','X')
      CALL NAME('','X')
      CALL NAME('Index Value','Y')    ! Axis name
      CALL plot_trees(r,num(p:q,1))  
      CALL tombox(1,cly(1),-1.2D0,1.2D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)  ! Gridlines
      CALL LINWID(1) ; CALL HEIGHT(22)
      CALL SETCLR(blue)
      CALL thickthin(r,wka(p:q),crn(p:q,1),num(p:q,1),7)
      CALL MESSAG("d) Chronology means 1.0",grl+50,grt+30)
      CALL SETCLR(red)
      CALL thickthin(r,wka(p:q),crn(p:q,11),num(p:q,11),7)
      CALL MESSAG("Chronology means reset",grl+900,grt+30)
      CALL SETCLR(black) ; CALL ENDGRF() 

      grt=1570 ; grb=1890
      p=1-cfy(1)+1 ; q=cyr(1) ; r=q-p+1 
      CALL TICKS(2,'Y')            ! No of ticks
      CALL LABELS('FLOAT','X')
      CALL LABDIG(-1,'Y')
      CALL NAME('Calendar Year','X')
      CALL NAME('Tree Count','Y')    ! Axis name
      CALL tombox(1,cly(1),0.D0,70.D0)
      CALL SETCLR(grey) ; CALL GRID(1,1)  ! Gridlines
      CALL LINWID(1) ; CALL HEIGHT(22)
      DO i=6,8
        CALL SETCLR(col(i-5))
        CALL thickthin(r,wka(p:q),DBLE(num(p:q,i+10)),num(p:q,i),7)
      ENDDO 
      CALL SETCLR(black)
      CALL MESSAG("e) Sub-sample chronology tree counts",grl+700,grt+28)
      CALL LABDIG(-1,'Y')
      CALL ENDGRF()
      CALL TICKS(5,'Y') ; CALL TICKS(10,'X')  ! No of ticks
      RETURN 
      END SUBROUTINE R2Fig5d
!-------------------------------------------------------------------
      SUBROUTINE R2Fig5() 
      IMPLICIT NONE                 
      INTEGER  :: i,j,k,p,q,r,w
      cnam(30)="R" 
      cnam(1)="../../raw/rcs/yml-all.raw" ; cnam(26)="Yamal Long"
      wnam(2)="Slowest Growth"
      wnam(3)="Medium Growth"
      wnam(4)="Fastest Growth"
      wnam(5)="All trees"
      wnam(6)="Slowest Growth"
      wnam(7)="Medium Growth"
      wnam(8)="Fastest Growth"
      CALL det_default() ; idb=2   ! Convert to normal
      CDsp=50 ; tst=4 ; sfo=2      ! Sort on growth rate 
      idt=-2 ; src=2 ; srcno=3     ! 3 curve RCS, SF=on  
      nc=0 ; CALL read_rft(cnam(1)) 
      DO k=1,11,10
        IF (k.NE.1) gtr=2        ! Reset means
        cf=k ; CALL detrend()  
        w=cyr(k) ; crn(1:w,k+1:k+7)=0.D0 ; num(1:w,k+1:k+7)=0
        CALL splinet(w,crn(1:w,k),CDsp,crn(1:w,k)) ! Smooth CRN
        DO j=2,5
          IF (j.EQ.5) THEN 
            i=mx
          ELSE
            i=j-1
          ENDIF
          p=sfy(i) ; q=sly(i) ; r=q-p+1
          cfy(j+k-1)=p ; cly(j+k-1)=q ; cyr(j+k-1)=r 
          crn(1:q,j+k-1)=mval(1:q,i)  ! Store RCS curves and counts
          num(1:q,j+k-1)=mcnt(1:q,i)
          okc(1:q,j+k-1)=mok(1:q,i)
          IF (j.LT.5) THEN
            crn(1:w,j+3+k)=xcsm(1:w,i)  ! Store Smoothed CRN and counts
            num(1:w,j+3+k)=xnum(1:w,i)
            okc(1:w,j)=cok(1:w,i)
          ENDIF
        ENDDO
      ENDDO 
      OPEN(19,FILE="R2fig5.col",IOSTAT=ios,STATUS="REPLACE")
      WRITE(19,'("Column  1 = Ring Age")')
      WRITE(19,'("Column  2 = Total sample count")')
      WRITE(19,'("Column  3 = Slowest Growth")')
      WRITE(19,'("Column  4 = Medium Growth")')
      WRITE(19,'("Column  5 = Fastest Growth")')
      WRITE(19,'("Column  6 = All trees")')
      DO i=1,cyr(5)
        WRITE(19,'(2I5,4F8.3)') i,num(i,5),crn(i,2:5)
      ENDDO

      WRITE(19,*)
      WRITE(19,'("Column  1 = Year")')
      WRITE(19,'("Column  2 = Total sample count")')
      WRITE(19,'("Column  3 = Slowest sample count")')
      WRITE(19,'("Column  4 = Medium sample count")')
      WRITE(19,'("Column  5 = Fastest sample count")')
      WRITE(19,'("Column  6 = Means 1.0 Slowest Growth")')
      WRITE(19,'("Column  7 = Means 1.0 Medium Growth")')
      WRITE(19,'("Column  8 = Means 1.0 Fastest Growth")')
      WRITE(19,'("Column  9 = Means reset Slowest Growth")')
      WRITE(19,'("Column 10 = Means reset Medium Growth")')
      WRITE(19,'("Column 11 = Means reset Fastest Growth")')
      WRITE(19,'("Column 12 = Chronology means 1.0")')
      WRITE(19,'("Column 13 = Chronology means reset")')
      DO i=1,cyr(1)
        WRITE(19,'(5I5,8F8.3)') cfy(1)-1+i,num(1,1),num(i,16:18), &
          crn(i,6:8),crn(i,16:18),crn(i,1),crn(i,11)
      ENDDO
      CLOSE(19)
      RETURN 
      END SUBROUTINE R2Fig5
!------------------------------------------------------------------------
     END MODULE rcs 

