;
; Plots low frequency MXD, TRW or ABD_MXD data versus temperature
; after calibration of the filtered series.
;
; Define periods for calibration, filters etc.
;
calper=[1881,1960]        ; calibration period
trv=0   ; 0=MXD 1=TRW 3=ABD_MXD (gets obs from MXD)
case trv of
  0: begin
    fnvar='mxd'
    end
  1: begin
    fnvar='trw'
    end
  3: begin
    fnvar='mxd'
    end
endcase
titadd=strupcase(fnvar)
;
; We want lots of overlapping band-passed filters, whose width increases
; linearly with time scale (from a width of 10 for the 0-10 year timescale,
; to a width of 50 for the 50-100 year timescale, noting that only an 80-yr
; run of data is analysed).
;
thalf=80
;
if trv eq 3 then begin
  ;
  ; Get all the age-banded MXD series (regions, hemi and back to 1400)
  ;
  restore,filename='/cru/u2/f055/treeharry/banding/bandregnh_fixed.idlsave'
  ; Gets: timey,nyr,nhts,allts,alleps,nreg,regname
  yrmxd=timey
  fullmxd=allts
  ; Combine hemispheric MXD with regional MXD
  nfull=n_elements(yrmxd)
  newfm=fltarr(nfull,nreg+2)*!values.f_nan
  newfm(*,0:nreg-1)=transpose(fullmxd(*,*))
  newfm(*,nreg)=nhts(*)
  newfm(*,nreg+1)=nhts(*)
  fullmxd=newfm
endif
;
; Get the hemispheric series
;
restore,filename='../compbest_'+fnvar+'_fixed1950.idlsave'
nhtemp=reform(comptemp(*,3))
alltemp=reform(comptemp(*,2))
nhmxd=reform(compmxd(*,2,0))
; Get rid of pre-1871 termperatures
knh=where(timey lt 1871)
nhtemp(knh)=!values.f_nan
alltemp(knh)=!values.f_nan
;
; Get regional series (temp, and truncated MXD)
;
restore,filename='../regbest_'+fnvar+'_fixed1950.idlsave'
; Gets: nreg,regname,regtemp,rawtemp,timey,bestmxd etc.
nyr=n_elements(timey)
;
; Add in the hemispheric series
;
regname=[regname,'ALL','NH']
newrt=fltarr(nyr,nreg+2)*!values.f_nan
newrt(*,0:nreg-1)=regtemp(*,*)
newrt(*,nreg)=alltemp(*)
newrt(*,nreg+1)=nhtemp(*)
regtemp=newrt
newrm=fltarr(nyr,nreg+2)*!values.f_nan
newrm(*,0:nreg-1)=bestmxd(*,*,1)
newrm(*,nreg)=nhmxd(*)
newrm(*,nreg+1)=nhmxd(*)
bestmxd=newrm
nreg=nreg+2
;
if trv eq 3 then begin
  ;
  ; Replace MXD with age-banded MXD (cutting out an appropriate period!)
  ;
  kmxd=where((yrmxd ge timey(0)) and (yrmxd le timey(nyr-1)),nmxd)
  bestmxd=fullmxd(kmxd,*)
  print,yrmxd(kmxd)
  ;
endif
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=4,ncol=3,layout='large'
if !d.name eq 'X' then begin
  wintim,ysize=850,xsize=700
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=14
endelse
def_1color,10,color='blue'
def_1color,11,color='red'
def_1color,12,color='mlgrey'
;
; Process each region separately
;
listreg=[0,1,2,3,4,6,7,8,9,10]
nlist=n_elements(listreg)
partit='('+['a','b','c','d','e','f','g','h','i','j','k']+') '
for jreg = 0 , nlist-1 do begin
  ireg=listreg(jreg)
  print,regname(ireg)
  ;
  ; Extract series
  ;
  tstemp=reform(regtemp(*,ireg))
  tsmxd=reform(bestmxd(*,ireg))
  ;
  ; Identify calibration and verification subsets
  ;
  calkl=where( finite(tstemp+tsmxd) and $
               (timey ge calper(0)) and (timey le calper(1)) , ncal )
  print,'Calibration period:'+string(calper,format='(2I5)')+$
    '  length:'+string(ncal,format='(I5)')
  ;
  ; Calibrate the MXD using a single band
  ;
  mkcalibrate,tsmxd(calkl),tstemp(calkl),fitcoeff=fitcoeff,autocoeff=autocoeff
  print,'RAW'
  print,'     r   alpha    beta SEalpha  SEbeta    sig    rmse'
  print,fitcoeff,format='(F6.2,4F8.4,F7.2,F8.4)'
  print,' a_mxd a_tem a_resid   ess'
  print,autocoeff,format='(3F6.2,F8.1)'
  sfac=fitcoeff(2)
  ;
  ; Calibrate the MXD, separately for each frequency band
  ;
  t1=reform(tstemp)
  m1=reform(tsmxd)
  filter_cru,/nan,thalf,tsin=t1,tslow=t2
  filter_cru,/nan,thalf,tsin=m1,tslow=m2
  mkcalibrate,m2(calkl),t2(calkl),fitcoeff=fitcoeff,autocoeff=autocoeff
  m4=m2(calkl)
  t4=t2(calkl)
  pred4=fitcoeff(1)+fitcoeff(2)*m4
  ;
  plot,timey(calkl),t4,$
    /xstyle,xtitle='Year',$
    ytitle='Temperature  (!Uo!NC)',$
    /ystyle,yrange=[-1,0.4],$
    title=partit(jreg)+regname(ireg)
  oplot,timey(calkl),pred4,thick=4
  oplot,!x.crange,[0,0],linestyle=1
  ;
endfor
;
end