;
; Reads in the age-banded ALL 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.
;
; *** Now averages bands 0-50 and 50-100, then averages all bands (with
; *** no weighting, but variance adjustment).
;
; *** Can optionally exclude any band with one 1 tree in it.
;
; *** Does it for fixed grids ***
;
fixedyr=1950
fnfix=string(fixedyr,format='(I4)')
;
allnyr=2000
allyr=findgen(allnyr)+1.
kper=where((allyr ge 1400) and (allyr le 1995),nyr)
timey=allyr(kper)
;
allts=fltarr(nyr)*!values.f_nan
rbarper=[1700,1994]
;
loadct,39
multi_plot,nrow=5
if !d.name eq 'X' then begin
  window,ysize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=12
endelse
;
; Read raw data
;
fn='mxd.allsites.BS5T1S1.17002000.sa.dat'+fnfix
;fn='/cru/keith1/harry/trees/neworder/mxd.allsites.all.BS5T1S1.17002000.pa.dat'
;fn='oldharryversion/mxd.allsites.all.BS5T1S1.17002000.sa.dat'
print,fn
;
; Check the number of columns
;
spawn,'head -1 '+fn+' | wc -w',osoutput
ncol=fix(osoutput)
ncol=ncol(0)
nband=(ncol-1)/3
;
if nband lt 5 then message,'Not enough bands with data!'
;
openr,1,fn
headlin=fltarr(ncol)
readf,1,headlin
rawdat=fltarr(ncol,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)*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
;ntree=rawdat(collist+2,*)
;misslist=where(ntree lt 2,nmiss)
;if nmiss gt 0 then tsdat(misslist)=!values.f_nan
; 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
tsnew=tsold
;
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
;
; Store in the 90 to 100 band to make following averaging simpler
; And make sure the ntree column of this band won't lead to removal of values.
;
rawdat(colts(9),*)=tsnew
rawdat(colts(9)+2,*)=100.
;
; Second stage of averaging: bands 0 to 50 yrs
;
print,'BANDS 0 TO 50 YEARS'
bandlist=indgen(5)+0
collist=colts(bandlist)
tsdat=rawdat(collist,*)
; Optionally remove those series made up from <2 trees
;ntree=rawdat(collist+2,*)
;misslist=where(ntree lt 2,nmiss)
;if nmiss gt 0 then tsdat(misslist)=!values.f_nan
; 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
tsnew=tsold
;
; Store in the 80 to 90 band to make following averaging simpler
; And make sure the ntree column of this band won't lead to removal of values.
;
rawdat(colts(8),*)=tsnew
rawdat(colts(8)+2,*)=100.
;
; Next stage of averaging: all bands
;
bandlist=indgen(nband-8)+8
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
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))
    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(float(allvar),1)
meanvar=totvar/float(numvar)
plot,timey,totvar,/xstyle,ytitle='Mean input variance'
;
print,'ALL BANDS'
tsdat=rawdat(collist,*)
; Optionally remove those series made up from <2 trees
;ntree=rawdat(collist+2,*)
;misslist=where(ntree lt 2,nmiss)
;if nmiss gt 0 then tsdat(misslist)=!values.f_nan
; 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)
tsnew=tsold
;
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
;
allts(*)=tsnew(*)
nhts=allts
;
y=allts
cpl_barts,timey,y,$
  /xstyle,xtitle='Year',$
  ytitle='Normalised MXD anomaly',$
  title='ALL SITES: '+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 save the data for subsequent analysis
;
save,filename='bandallnh_science_fixed'+fnfix+'.idlsave',$
  timey,nyr,nhts
;
end