;
; Calibrates and plots MXD against temperature
; DOES IT FOR AGE-BANDED DENSITY FROM ALLBANDS
;
; Define periods for calibration
;
calper=[1881,1960]        ; calibration period
verper=[1961,1994]        ; verification period
thalf=10.
donh=1                    ; 0=do NH 1=do ALL
;
if donh eq 0 then nhtit='NH' else nhtit='ALL'
;
; Get the hemispheric temperature series
;
restore,filename='compbest_mxd_fixed1950.idlsave'
if donh eq 0 then nhtemp=reform(comptemp(*,3)) $
             else nhtemp=reform(comptemp(*,2))
; Get rid of pre-1871 temperatures
knh=where(timey lt 1871)
nhtemp(knh)=!values.f_nan
nyrtemp=n_elements(timey)
yrtemp=timey
;
; Get all the age-banded MXD series (regions, hemi and back to 1400)
;
restore,filename='/cru/u2/f055/treeharry/banding/bandallnh_fixed1950.idlsave'
;restore,filename='/cru/u2/f055/treeharry/banding/bandallnh_science_fixed1950.idlsave'
; Gets: timey,nyr,nhts
yrmxd=timey
nhmxd=nhts
nyrmxd=nyr
;
; Identify overlap between instrumental and MXD series
;
kcomp=where( (yrmxd ge yrtemp(0)) and (yrmxd le yrtemp(nyrtemp-1)) )
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=3,layout='large',/landscape
if !d.name eq 'X' then begin
  window,xsize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=9
endelse
def_1color,10,color='blue'
def_1color,11,color='red'
;
range0=-1.5
range1=1.0
;
calregts=fltarr(nyrmxd)*!values.f_nan
;
  ; Extract series
  ;
  tstemp=nhtemp
  tsmxd=nhmxd(kcomp)
  tsfull=nhmxd
  ;
  ; Identify calibration and verification subsets
  ;
  calkl=where( finite(tstemp+tsmxd) and $
               (yrtemp ge calper(0)) and (yrtemp le calper(1)) , ncal )
  print,'Calibration period:'+string(calper,format='(2I5)')+$
    '  length:'+string(ncal,format='(I5)')
  ;
  ; Calibrate the MXD
  ;
dummy=where(finite(tsmxd),nnnnnn)
if nnnnnn gt 2 then begin
  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)'
  ;
  ; And low and high passed too!
  ;
  filter_cru,/nan,thalf,tsin=tsmxd,tshigh=mxdhi,tslow=mxdlo
  filter_cru,/nan,thalf,tsin=tstemp,tshigh=temphi,tslow=templo
  mkcalibrate,mxdhi(calkl),temphi(calkl),fitcoeff=hicoeff,autocoeff=autohi
  print,'HIGH-PASSED',thalf
  print,'     r   alpha    beta SEalpha  SEbeta    sig    rmse'
  print,hicoeff,format='(F6.2,4F8.4,F7.2,F8.4)'
  print,' a_mxd a_tem a_resid   ess'
  print,autohi,format='(3F6.2,F8.1)'
  mkcalibrate,mxdlo(calkl),templo(calkl),fitcoeff=locoeff,autocoeff=autolo,$
    nlag=5
  print,'LOW-PASSED',thalf
  print,'     r   alpha    beta SEalpha  SEbeta    sig    rmse'
  print,locoeff,format='(F6.2,4F8.4,F7.2,F8.4)'
  print,' a_mxd a_tem a_resid   ess'
  print,autolo,format='(3F6.2,F8.1)'
  ;
  ; Generate calibrated record and uncertainties etc.
  ; Use the extended MXD record (after checking it matches the short one!)
  ;
  xxx=tsfull(kcomp)
  yyy=tsmxd
  tse=total( (xxx-yyy)^2 , /nan )
;  print,tse
  if tse gt 0.03 then message,'Series do not match'
  calmxd=tsfull*fitcoeff(2)+fitcoeff(1)
  calregts(*)=calmxd(*)
  filter_cru,/nan,thalf,tsin=tsfull,tshigh=fullhi,tslow=fulllo
  sepmxd=(fullhi*hicoeff(2)+hicoeff(1))+(fulllo*locoeff(2)+locoeff(1))
  ;
  ; For overlap period with non-missing data, compute correlation between
  ; the temperature series and the 2-band reconstruction
  ;
  kkk=calkl+(yrtemp(0)-yrmxd(0))
  print,'r(Two-band calibrated series,Temp)='+$
    string(correlate(tstemp(calkl),sepmxd(kkk)),format='(F6.2)')
  ;
  plot,yrtemp,tstemp,$
    /xstyle,xrange=[1400,1994],xtitle='Year',$
    ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
    /ystyle,yrange=[range0,range1],$
    title=nhtit
  oplot,yrmxd,calmxd,color=10       ;,thick=3
  oplot,yrmxd,sepmxd,color=11       ;,thick=3
  oplot,!x.crange,[0.,0.],linestyle=1
  filter_cru,/nan,thalf,tsin=tstemp,tslow=ylow
  oplot,yrtemp,ylow,thick=3
  filter_cru,/nan,thalf,tsin=calmxd,tslow=ylow
  oplot,yrmxd,ylow,thick=3,color=10
  ;
  plot,yrtemp,tstemp,thick=2,$
    /xstyle,xrange=[1850,2000],xtitle='Year',$
    ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
    title=nhtit
  oplot,yrmxd,calmxd,color=10,thick=3
  oplot,yrmxd,sepmxd,color=11,thick=3
  oplot,!x.crange,[0.,0.],linestyle=1
endif
;
; Now output the calibrated series
;
nyr=nyrmxd
timey=yrmxd
nhtemp=calregts
save,filename='bandalltemp_calibrated.idlsave',$
  nyr,nhtit,timey,nhtemp
;
end