;
; From calibrated regional MXD series, computes residual from the
; temperature series, and correlates this against pre-computed
; regional TOMS ozone series for 1978-1993.
;
doclassic=1     ; 0=inverse approach, 1=classical approach
;
doplot=[1,0,1,0,1,0,0,0,0]
;doplot=intarr(9)+1
;
; Get the calibrated MXD and temperature series
;
restore,filename='regtemp_calibrated.allversion.idlsave'
; Gets:  nyr,nreg,calregts,regname,timey,tempregts,tempnyr,temptimey
;
; Extract the instrumental period (and remove the ALL/NH region)
;
kmxd=where(timey ge temptimey(0),nmxd)
ktem=where(temptimey le timey(nyr-1),ntem)
if nmxd ne ntem then message,'Oooops!'
calregts=calregts(kmxd,0:nreg-2)
tempregts=tempregts(ktem,0:nreg-2)
temptimey=temptimey(ktem)
;
; NEUR regional mean in 1993 is based on just one site, so set to missing!
;
gotreg=where(regname eq 'NEUR',ngotr)
gotyr=where(temptimey eq 1993,ngoty)
if (ngotr eq 1) and (ngoty eq 1) then begin
  print,'DROPPING 1993 FROM NEUR'
  calregts[gotyr[0],gotreg[0]]=!values.f_nan
endif
;
; Now get the ozone series
;
restore,filename='/cru/u2/f055/data/sat/toms/monthlyozone/reg_ozone.idlsave'
; Gets:  nreg,regname,regts,timey,nyr,nmon,monname
;
; Find overlap period
;
kover=where( (temptimey ge timey(0)) and (temptimey le timey(nyr-1)) , nover)
;
; Now prepare for plotting
;
loadct,39
multi_plot,nrow=6
if !d.name eq 'X' then begin
  window,ysize=800
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=12
endelse
def_1color,20,color='red'
def_1color,21,color='green'
;
for ireg = 0 , nreg-1 do begin
 if doplot(ireg) eq 1 then begin
  ;
  ; If we're adopting the "classical" regression approach, we need to reverse
  ; the direction of the regression, and no attempt to predict the MXD series
  ; from the temperatures, and subtract this prediction from the MXD series
  ; to leave the residual MXD series!
  ;
  if doclassic ne 0 then begin
    xxmxd=reform(calregts[*,ireg])
    xxtem=reform(tempregts[*,ireg])
    xxyr=temptimey
    mknormal,xxmxd,xxyr,refperiod=[1881,1960]
    calregts[*,ireg]=xxmxd[*]
    xxkl=where((xxyr ge 1881) and (xxyr le 1960) and $
      (finite(xxmxd+xxtem)),nkeep)
    print,regname[ireg]
    print,nkeep
    xxmxd=xxmxd[xxkl]
    xxtem=xxtem[xxkl]
    rdt=correlate(xxmxd,xxtem)
    print,rdt
    xxcoeff=linfit(xxtem,xxmxd)
    print,xxcoeff
    tempregts[*,ireg]=xxcoeff[0]+xxcoeff[1]*tempregts[*,ireg]
  endif
  ;
  if doclassic eq 0 then begin
    yran=[-2.7,2.1]
    ytit='Temperature (!Uo!NC wrt 1961-90)'
    yran2=[-1.8,1.6]
    ytit2='Temperature residual (!Uo!NC)'
  endif else begin
    yran=[-3.5,3.5]
    ytit='Normalised tree-ring density'
    yran2=[-3.5,3.5]
    ytit2='Density residual'
  endelse
  ;
  pause
  plot,temptimey,calregts(*,ireg),thick=3,$
    /xstyle,xrange=[1860,2000],$
    /ystyle,yrange=yran,ytitle=ytit,$
    ymargin=[0,2],title=regname(ireg),xtickformat='nolabels'
  oplot,temptimey,tempregts(*,ireg),thick=3,color=20
  oplot,!x.crange,[0.,0.],linestyle=1
  if doclassic ne 0 then begin
    xyouts,align=1,1995,yran[0]+0.1*(yran[1]-yran[0]),$
      'r ='+string(rdt,format='(2F5.2)')
  endif
  ;
  x1a=calregts(*,ireg)
  x1b=tempregts(*,ireg)
  if doclassic eq 0 then begin
    x1=x1a-x1b
  endif else begin
    x1=x1a-x1b
  endelse
  x2=reform(regts(ireg,3,*))
;  x2=total(reform(regts(ireg,2:3,*)),1)/2.
  mknormal,x2
  x2=x2*0.6
  plot,temptimey,x1,thick=3,$
    /xstyle,xrange=[1860,2000],xtitle='Year',$
    /ystyle,yrange=yran2,ytitle=ytit2,$
    ymargin=[4,0]
; Arbitrary lag (no physical basis!)
;  nx2=n_elements(x2)
;  x2=shift(x2,-1) & x2[nx2-1]=!values.f_nan
  oplot,timey,x2,thick=3,color=21
  oplot,!x.crange,[0.,0.],linestyle=1
  r3=mkcorrelation(x1(kover),x2,datathresh=5,filter=5.)
  print,'T vs O3'
  print,mkcorrelation(x1b(kover),x2,datathresh=5,filter=5.)
  print,'MXD vs O3'
  print,mkcorrelation(x1a(kover),x2,datathresh=5,filter=5.)
  print,'MXDresidual vs O3'
  print,mkcorrelation(x1(kover),x2,datathresh=5,filter=5.)
  ;
;  if finite(r3(0)) then xyouts,1960,1,'r ='+string(r3[0:1],format='(2F5.2)')
  if finite(r3(0)) then begin
    xyouts,align=1,1995,yran[0]+0.1*(yran[1]-yran[0]),$
      'r ='+string(r3[0],format='(2F5.2)')
  endif
  ;
 endif
endfor
;
end