;
; Reads in the age-banded regional series produced by Harry's program
; saband10.f.  For each region, averages across bands 50yr-100yr, then
; across bands 50yr-300yr.  Both these averaging processes are variance-
; adjusted, and from the second one the rbar is used to compute a time-
; dependent EPS.  The regional mean series are then averaged to form a
; northern hemisphere series, weighted by their EPS, and variance adjusted.
;
dodata=1      ; 0=do not, 1=do output data to ascii file
;
; EPS weighting can optionally be turned off
epsw=1        ; 0=unweighted, 1=weighted by EPS
;
; *** Can optionally exclude any band with one 1 tree in it.
mintree=2     ; 1 or <1 means include all, >1 means include only bands with
              ; >=mintree in them
;
; *** Does it for fixed grids ***
fixedyr=1950        ; 1400, 1500, 1600, 1700 ,1800, 1950, 1988
fnfix=string(fixedyr,format='(I4)')
if fixedyr eq 1950 then fnfix=''
;
; Define components of the output file names
;
fnband='_50-300'
fnband='_20-550'
if mintree gt 1 then fnband=fnband+'_min'+string(mintree,format='(I1)')
if epsw ne 0 then fnband=fnband+'_eps'
;
if dodata then begin
missval=-999.999
fnout='bandregnh'+fnband+'_fixed'+fnfix+'.dat'
print,fnout
openw,2,fnout
;
printf,2,'Data from Tim Osborn, April 2004'
printf,2
printf,2,'Age-band decomposition (ABD) results applied to the'
printf,2,'Schweingruber tree-ring density data set.'
printf,2
printf,2,'Only trees in age bands between 50 and 300 years old are used.'
if mintree gt 1 then printf,2,'Only age bands with >=',mintree,' tree cores in them are used.'
printf,2
printf,2,'Time series are given for each of 9 regions, showing the'
printf,2,'the number of age bands with data in them, the (rbar) average'
printf,2,'inter-band correlation, the EPS, and the'
printf,2,'regional-average tree-ring density.'
endif
;
; Define region names.  Names are different to those I use, but order
; should be the same to make later comparison easier.
;
regname=['northwestern.europe',$
         'southwestern.europe',$
         'northern.siberia',$
         'eastern.siberia',$
         'centralsiberia',$
         'tibetianplateau',$
         'western.northamerica',$
         'northwestern.canada',$
         'centraleastern.canada']
nreg=n_elements(regname)
;
myreg=['NEUR','SEUR','NSIB','ESIB','CAS','TIBP','WNA','NWNA','ECCA']
;
; Specify the number of bands used in each region
;
; **These values are now overwritten, because they're lower for the early
; **fixed grids!
nband=[19,20,21,24,18,20,22,17,17]
ncol=nband*3 + 1
;
allnyr=2000
allyr=findgen(allnyr)+1.
kper=where((allyr ge 1400) and (allyr le 1995),nyr)
timey=allyr(kper)
;
allts=fltarr(nreg,nyr)*!values.f_nan
alleps=fltarr(nreg,nyr)*!values.f_nan
rbarper=[1700,1994]
;
loadct,39
multi_plot,nrow=5
if !d.name eq 'X' then begin
  window,ysize=750
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=12
endelse
;
; Repeat for each region
;
for ireg = 0 , nreg-1 do begin
  ;
  ; Read raw data
  ;
  fn='mxd.'+regname(ireg)+'.all.BS5T1S1.17002000.sa.dat'+fnfix
  print,fn
  ;
  ; Check the number of columns
  ;
  spawn,'head -1 '+fn+' | wc -w',osoutput
  ncol(ireg)=fix(osoutput)
  nband(ireg)=(ncol(ireg)-1)/3
  ;
  if nband(ireg) lt 5 then begin
    print,'Skipping region!'
  endif else begin
  ;
  openr,1,fn
  headlin=fltarr(ncol(ireg))
  readf,1,headlin
  rawdat=fltarr(ncol(ireg),allnyr)
  readf,1,rawdat
  close,1
  ;
  ; Just keep the 1400 to 1994 portions
  ;
  rawdat=rawdat(*,kper)
  ;
  ; Locate missing data
  ;
  ml=where(rawdat eq 9990.,nmiss)
  rawdat(ml)=!values.f_nan
  ;
  colts=indgen(nband(ireg))*3+1
  ;
  ; First stage of averaging: bands 50 to 100 yrs
  ;
  print,'BANDS 50 TO 100 YEARS'
  bandlist=indgen(5)+5             ; was (5)+5
  collist=colts(bandlist)
  tsdat=rawdat(collist,*)
  ; Optionally remove those series made up from <2 trees
  if mintree gt 1 then begin
    ntree=rawdat(collist+2,*)
    misslist=where(ntree lt mintree,nmiss)
    if nmiss gt 0 then tsdat(misslist)=!values.f_nan
  endif
  ; End of option
  totdat=total(tsdat,1,/nan)
  numdat=total(finite(tsdat),1)
  tsold=totdat/float(numdat)
  tsnew=var_adjust(tsdat,timey,refperiod=rbarper,rbar=fullrbar,/mkrbar,$
    /normalise)
  ; If 1700-1994 had insufficient data to compute rbar, use full period
  if finite(fullrbar) eq 0 then begin
    tsnew=var_adjust(tsdat,timey,rbar=fullrbar,/mkrbar,/normalise)
  endif
  ;
  pause
  plot,timey,numdat,/xstyle,/ystyle,yrange=[0,6],$
    title='rbar='+string(fullrbar,format='(F6.2)')
  neff=float(numdat)/(1.+(float(numdat)-1.)*fullrbar)
  oplot,timey,neff,thick=4
  legend,['Number of 50-100 decadal bands','Effective independent bands'],$
    thick=[1,4],charsize=0.5
  ;
  ; Store in the 90 to 100 band to make following averaging simpler
  ; And make sure the ntree column for this band won't lead to the removal of
  ; any data.
  ;
  rawdat(colts(9),*)=tsnew
  rawdat(colts(9)+2,*)=100.
  ;
  ; Next stage of averaging: bands 50 to 300 yrs
  ;
  bandlist=indgen(5)+9              ; was (5)+9
  collist=colts(bandlist)
  nuse=n_elements(collist)
  ;
  ; Compute and plot running variance of 5 input series
  ;
  print,'COMPUTING INPUT VARIANCE'
  allvar=fltarr(nuse,nyr)*!values.f_nan
  wlen=50
  ; We normalise the input series prior to averaging them, so lets see what
  ; the running variances are like *after* normalising them
  normdat=rawdat[collist,*]
  mknormal,normdat,timey,refperiod=rbarper
  for iyr = 0 , nyr-wlen do begin
    print,timey(iyr),format='($,I4)'
    for iband = 0 , nuse-1 do begin
;      oneseg=reform(rawdat(collist(iband),iyr:iyr+wlen-1))
      oneseg=reform(normdat(iband,iyr:iyr+wlen-1))
      kl=where(finite(oneseg),nkeep)
      if nkeep ge 10 then begin
        print,'.',format='($,A1)'
        onemon=moment(oneseg(kl))
        allvar(iband,iyr)=onemon(1)
      endif
    endfor
  endfor
  print
  totvar=total(allvar,1,/nan)
  numvar=total(finite(allvar),1)
  meanvar=totvar/float(numvar)
  plot,timey,meanvar,/xstyle,ytitle='Variance',thick=4,$
    /ystyle,yrange=[0,max(allvar,/nan)],$
   title='Variance of the input band time series'
  for iband = 0 , nuse-1 do oplot,timey,allvar[iband,*]
  ;
  print,'BANDS 50 TO 300 YEARS'
  tsdat=rawdat(collist,*)
  ; Optionally remove those series made up from <2 trees
  if mintree gt 1 then begin
    ntree=rawdat(collist+2,*)
    misslist=where(ntree lt mintree,nmiss)
    if nmiss gt 0 then tsdat(misslist)=!values.f_nan
  endif
  ; End of option
  totdat=total(tsdat,1,/nan)
  numdat=total(finite(tsdat),1)
  tsold=totdat/float(numdat)
  tsnew=var_adjust(tsdat,timey,refperiod=rbarper,rbar=fullrbar,/mkrbar,$
    /normalise)
  ;
  plot,timey,numdat,/xstyle,/ystyle,yrange=[0,6],$
    title='rbar='+string(fullrbar,format='(F6.2)')
  neff=float(numdat)/(1.+(float(numdat)-1.)*fullrbar)
  oneeps=fullrbar / (fullrbar + (1.-fullrbar)/float(numdat) )
  oplot,timey,neff,thick=4
  legend,['Number of 50-300 decadal bands','Effective independent bands'],$
    thick=[1,4],charsize=0.5
  ;
  plot,timey,oneeps,/xstyle,/ystyle,yrange=[0.,1.],$
    title='EPS'
  legend,['EPS based on n & overall rbar','Time-varying rbar',$
    'EPS based on n & time-varying rbar'],$
    thick=[1,3,5],charsize=0.5,color=[!p.color,!p.color,240]
  ;
  ; Also compute a time-dependent correlation matrix, and then rbar(t)
  ; and then a new eps(t).  If rbar becomes very low, or if it becomes
  ; zero simply because there is only one band with data, then we use
  ; the oldest positive rbar value.
  ;
  tdrbar=fltarr(nyr)
  print,'Century rbar:'
  for icent = 0 , 11 do begin
    stcent=float(icent)*50.+1400.
    print,stcent,format='($,I5)'
    kcent=where((timey ge stcent) and (timey le stcent+49.))
    centts=tsdat(*,kcent)
    tdrbar(kcent)=mkrbar(centts) > 0.
  endfor
  for iyr = nyr-100 , 0 , -1 do begin
    if numdat(iyr) le 0. then begin
      tdrbar(iyr)=!values.f_nan
    endif else begin
      if finite(tdrbar(iyr)) eq 0 then tdrbar(iyr)=tdrbar(iyr+1) $
      else if tdrbar(iyr) le 0. then tdrbar(iyr)=tdrbar(iyr+1)
    endelse
  endfor
  tdeps=tdrbar / (tdrbar + (1.-tdrbar)/float(numdat) )
  ;
  ; Filter (optional) the EPS series so that it doesn't suddenly change
  ; just because a tree moves from one age band to another
  ;
  ; Only worth doing if we're actually going to use the EPS for later weighting
  if epsw ne 0 then begin
    print,'Smoothing EPS'
    filter_cru,20.,/nan,tsin=tdeps,tslow=tdepslow
    tdeps=tdepslow
  endif
  ;
  oplot,timey,tdeps,thick=3,color=240
  oplot,timey,tdrbar,thick=5
  ;
  ; Now output the data
  ;
  if dodata then begin
  tdrbaro=tdrbar
  tdepso=tdeps
  tsnewo=tsnew
  ml=where(finite(tdrbaro) eq 0,nmiss)
  if nmiss gt 0 then tdrbaro[ml]=missval
  ml=where(finite(tdepso) eq 0,nmiss)
  if nmiss gt 0 then tdepso[ml]=missval
  ml=where(finite(tsnewo) eq 0,nmiss)
  if nmiss gt 0 then tsnewo[ml]=missval
  printf,2
  printf,2,myreg[ireg],' / ',regname[ireg]
  printf,2
  printf,2,'Year','Nband','Rbar','EPS','Tree-density',format='(5A15)'
  for iyr = 0 , nyr-1 do begin
    printf,2,timey[iyr],numdat[iyr],tdrbaro[iyr],tdepso[iyr],tsnewo[iyr],format='(I15,4F15.3)'
  endfor
  endif
  ;
; OPTION TO GIVE UNIFORM WEIGHTING OF ALL REGIONS!!!!!
  if epsw eq 0 then tdeps(*)=1.
  ;
  allts(ireg,*)=tsnew(*)
  alleps(ireg,*)=tdeps(*)
  ;
  y=reform(allts(ireg,*))
  cpl_barts,timey,y,$
    /xstyle,xtitle='Year',$
    ytitle='Normalised MXD anomaly',$
    title=regname(ireg)+fnfix
  filter_cru,/nan,20.,tsin=y,tslow=ylow
  oplot,timey,ylow,thick=4
  oplot,!x.crange,[0.,0.]
  ;
  endelse
  ;
endfor
;
; Average across regions to give NH
;
nhtot=total(allts,1,/nan)
nhnum=total(finite(allts),1)
nhold=nhtot/float(nhnum)
;
; Compute rbar between the regions for each 50-yr period, to measure any
; deteriation.
;
rbarvary=fltarr(nyr)
print,'Changing rbar:'
for icent = 0 , 11 do begin
  stcent=float(icent)*50.+1400.
  print,stcent,format='($,I5)'
  kcent=where((timey ge stcent) and (timey le stcent+49.))
  centts=allts(*,kcent)
  rbarvary(kcent)=mkrbar(centts)
endfor
for iyr = nyr-100 , 0 , -1 do begin
  if nhnum(iyr) le 0. then begin
    rbarvary(iyr)=!values.f_nan
  endif else begin
    if finite(rbarvary(iyr)) eq 0 then rbarvary(iyr)=rbarvary(iyr+1)
  endelse
endfor
;
nhnew=var_adjust(allts,timey,weight=alleps,$
  refperiod=rbarper,rbar=fullrbar,/mkrbar,/normalise,/td_rbar)
;
nhts=nhnew
;
; Plot it
;
pause
;plot,timey,rbarvary,thick=3,title='Varying rbar between regions',$
;  yrange=[-0.2,0.5],/xstyle,/ystyle
;oplot,!x.crange,[0.,0.],linestyle=1
;
;plot,timey,timey,/nodata,/xstyle,/ystyle,yrange=[0.,1.],$
;  title='Time-dependent EPS'
;for ireg = 0 , nreg-1 do oplot,timey,alleps(ireg,*),thick=ireg+1
;
plot,timey,fullrbar,/xstyle,/ystyle,yrange=[0.,0.5],$
  title='Time-dependent rbar (due to regions dropping out, but not due to deterioration back in time)'
;
plot,timey,nhnum,/xstyle,/ystyle,yrange=[0,10]
neff=float(nhnum)/(1.+(float(nhnum)-1.)*fullrbar)
oplot,timey,neff,thick=4
legend,['Number of regions with data','Effective independent regions'],$
  thick=[1,4],charsize=0.5
;
; Compute EPS, but note that early rbar is set to 1 when only one region
for i = 10 , 0 , -1 do begin
  if fullrbar[i] eq 1. then fullrbar[i]=fullrbar[i+1]
endfor
oneeps=fullrbar / (fullrbar + (1.-fullrbar)/float(nhnum) )
plot,timey,oneeps,/xstyle,/ystyle,yrange=[0,1],$
  title='EPS'
;
y=nhts
cpl_barts,timey,y,$
  /xstyle,xtitle='Year',$
  ytitle='Normalised MXD anomaly',$
  title='NH mean using age-banded MXD, fixed grid: '+fnfix
filter_cru,/nan,20.,tsin=y,tslow=ylow
oplot,timey,ylow,thick=4
oplot,!x.crange,[0.,0.]
;
print,'COMPUTING OUTPUT VARIANCE'
nhvar=fltarr(nyr)*!values.f_nan
wlen=50
for iyr = 0 , nyr-wlen do begin
  print,timey(iyr),format='($,I4)'
  oneseg=nhts(iyr:iyr+wlen-1)
  kl=where(finite(oneseg),nkeep)
  if nkeep ge 10 then begin
    onemon=moment(oneseg(kl))
    nhvar(iyr)=onemon(1)
  endif
endfor
print
;
plot,timey,nhvar,/xstyle,ytitle='NH variance',$
  xtitle='Start year of 50-yr sliding window'
;
; Now output the data
;
if dodata then begin
tdrbaro=fullrbar
tdepso=oneeps
tsnewo=nhts
ml=where(finite(tdrbaro) eq 0,nmiss)
if nmiss gt 0 then tdrbaro[ml]=missval
ml=where(finite(tdepso) eq 0,nmiss)
if nmiss gt 0 then tdepso[ml]=missval
ml=where(finite(tsnewo) eq 0,nmiss)
if nmiss gt 0 then tsnewo[ml]=missval
printf,2
printf,2,'Average across all regions (i.e., a quasi-NH series)'
printf,2
printf,2,'Year','Nregion','Rbar','EPS','Tree-density',format='(5A15)'
for iyr = 0 , nyr-1 do begin
  printf,2,timey[iyr],nhnum[iyr],tdrbaro[iyr],tdepso[iyr],tsnewo[iyr],format='(I15,4F15.3)'
endfor
endif
;
close,2
;
; Now save the data for subsequent analysis
;
save,filename='bandregnh'+fnband+'_fixed'+fnfix+'.idlsave',$
  timey,nyr,nhts,allts,alleps,nreg,regname
;
end