;
; 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.
;
; *** Now averages bands 20-100yr first (not just 50-100yr), although they
; *** still only count as a single band.  And also uses bands 300-550yr,
; *** after first averaging them into a single band (weighting & adjusted etc.)
;
; *** Can optionally exclude any band with one 1 tree in it.
;
; *** Does it for fixed grids ***
;
fixedyr=1950
fnfix=string(fixedyr,format='(I4)')
if fixedyr eq 1950 then fnfix=''
;
; 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)
;
; 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=850
  !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(8)+2             ; 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
    print,'Re-doing variance adjustment using full period to compute rbar'
    tsnew=var_adjust(tsdat,timey,rbar=fullrbar,/mkrbar,/normalise)
  endif
  ;
  pause
  plot,timey,numdat,/xstyle,/ystyle,yrange=[0,8],$
    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
  ; Also replace the ntree column of 90-100 to make sure that no data is
  ; subsequently removed!
  ;
  rawdat(colts(9),*)=tsnew
  rawdat(colts(9)+2,*)=100.
  ;
  ; Second stage of averaging: bands 300 to 550 yrs
  ;
  print,'BANDS 300 TO 550 YEARS'
  bandlist=indgen(5)+14
  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
    print,'Re-doing variance adjustment using full period to compute rbar'
    tsnew=var_adjust(tsdat,timey,rbar=fullrbar,/mkrbar,/normalise)
  endif
  ;
  ; Store in the 300 to 350 band to make following averaging simpler
  ;
  rawdat(colts(14),*)=tsnew
  rawdat(colts(14)+2,*)=100.
  ;
  ; Next stage of averaging: bands 50 to 300 yrs
  ;
  bandlist=indgen(6)+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
  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 25 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,'BANDS 50 TO 300 YEARS'
  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)
  ;
  plot,timey,numdat,/xstyle,/ystyle,yrange=[0,7],$
    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
  ;
  plot,timey,oneeps,/xstyle,/ystyle,yrange=[0.,1.]
    title='EPS'
  ;
  ; 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
  ;
  print,'Smoothing EPS'
  filter_cru,20.,/nan,tsin=tdeps,tslow=tdepslow
  tdeps=tdepslow
; OPTION TO GIVE UNIFORM WEIGHTING OF ALL REGIONS!!!!!
; tdeps(*)=1.
  ;
  oplot,timey,tdeps,thick=3
  oplot,timey,tdrbar,thick=5
  ;
  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)'
;
plot,timey,nhnum,/xstyle,/ystyle,yrange=[0,10]
neff=float(nhnum)/(1.+(float(nhnum)-1.)*fullrbar)
oplot,timey,neff,thick=4
;
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 save the data for subsequent analysis
;
save,filename='bandregnh_fixed'+fnfix+'.idlsave',$
  timey,nyr,nhts,allts,alleps,nreg,regname
;
end