;
; Plots warm & cold season temperatures from paleo, instr and model
;
thalf=50.       ; filter for high or low-pass
;
loadct,39
multi_plot,nrow=3,layout='large'
if !d.name eq 'X' then begin
  wintim,ysize=900
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,bold=1,font_size=20
endelse
def_1color,20,color='red'
def_1color,21,color='deepblue'
;
; Get HadCM2
;
;restore,'/cru/u2/f055/nao/ddc/science_hc2_warm&cold.idlsave'
; Gets: timey,ts1,ts2,ts3
;modyr=timey
;modcold=ts1
;modwarm=ts2
restore,'/cru/u2/f055/data/mod/hc_hadcm2/sulensemblemean_tave.idlsave'
; Gets: timey,nyr,tavets[4,nyr]
modyr=timey
modcold=reform(tavets[3,*])
modwarm=reform(tavets[2,*])
;
; Get observations
;
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/crutem2_18512001.mon.nc')
tsmon=crunc_rddata(ncid,tst=[1851,0],tend=[2001,11],grid=gtemp)
ncdf_close,ncid
nmon=12
ntemp=gtemp.nt
nyrtemp=ntemp/nmon
obsyr=reform(gtemp.year,nmon,nyrtemp)
obsyr=reform(obsyr(0,*))
; Compute the northern hemisphere >20N land series
; First extract >20N rows
kl=where(gtemp.y gt 20.)
ylat=gtemp.y(kl)
tsnorth=tsmon(*,kl,*)
; Compute latitude-weighted mean
nhmon=globalmean(tsnorth,ylat)
; Compute seasonal/annual mean
nhmon=reform(nhmon,nmon,nyrtemp)
obswarm=mkseason(nhmon,3,8,datathresh=3)
obscold=mkseason(nhmon,9,2,datathresh=3)
;
; Read Esper's raw reconstruction, then scale to either warm or cold season
openr,1,'/cru/u2/f055/data/paleo/esper2002/esper.txt'
readf,1,nyr
headst=''
readf,1,headst
rawdat=fltarr(7,nyr)
readf,1,rawdat
close,1
palyr=reform(rawdat(0,*))
ts=reform(rawdat(1,*))
; Convert from raw values to deg C relative to 1961-1990
palwarm=ts*1.7138-2.1161       ; via regression vs. NH land>20 warm
palcold=ts*2.1361-2.7758       ; via regression vs. NH land>20 cold
;
; Now plot the results
;
filter_cru,thalf,/nan,tsin=modwarm,tslow=modwarml
filter_cru,thalf,/nan,tsin=modcold,tslow=modcoldl
filter_cru,thalf,/nan,tsin=obswarm,tslow=obswarml
filter_cru,thalf,/nan,tsin=obscold,tslow=obscoldl
filter_cru,thalf,/nan,tsin=palwarm,tslow=palwarml
filter_cru,thalf,/nan,tsin=palcold,tslow=palcoldl
  ;
plot,palyr,palwarml,/nodata,$
  xstyle=1,xrange=[830,2100],xtitle='Year  (AD)',$
  /ystyle,ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
  yrange=[-1.2,0.5]
oplot,!x.crange,[0.,0.],linestyle=1
  ;
oplot,modyr,modwarml,thick=2,color=20
oplot,modyr,modcoldl,thick=2,color=21
oplot,obsyr,obswarml,thick=6,color=20
oplot,obsyr,obscoldl,thick=6,color=21
oplot,palyr,palwarml,thick=2,color=20
oplot,palyr,palcoldl,thick=2,color=21
;
oplot,[!x.crange[0],1990],[2.2,2.2]
xyouts,0.5*(!x.crange[0]+1990),1.75,align=0.5,'palaeo',charsize=0.8
oplot,[1851,2001],[3.2,3.2]
xyouts,0.5*(1851+2001),2.75,align=0.5,'instr',charsize=0.8
oplot,[1860,!x.crange[1]],[4.2,4.2]
xyouts,0.5*(!x.crange[1]+1860),3.75,align=0.5,'model',charsize=0.8
  ;
end