;
; Computes correlation maps between JJA gridded temperatures and the
; northern hemisphere mean JJA temperature and the residual of the
; NH JJA temperature minus Phil's NH10 reconstruction.
;
dojja=0     ; 0=Apr-Sep  1=Jun-Aug
if dojja eq 0 then titadd='Apr-Sep' else titadd='Jun-Aug'
;
; Read in Phil's series and adjust appropriately to degrees C wrt 61-90
;
print,'Reading in NH10 series'
perst=1881
peren=1980
openr,1,'phil_nhrecon.dat'
nyr=992
rawdat=fltarr(4,nyr)
readf,1,rawdat,format='(I5,F7.2,I3,F7.2)'
close,1
timey=reform(rawdat(0,*))
ts=reform(rawdat(3,*))
kl=where((timey ge perst) and (timey le peren),nyr)
timey=timey(kl)
ts=ts(kl)
; Convert from normalised values to deg C relative to 1961-1990
ts=ts*0.521 - 0.1134
;
nh10yr=timey
nh10ts=ts
;
; Now read in the gridded land and marine dataset
;
nmon=12
yrstart=1881
yrend=1980
nyr=yrend-yrstart+1
fn1='/cru/u2/f055/data/obs/grid/surface/ist_ipcc_18561999.mon.nc'
timey=findgen(nyr)+yrstart
;
; Read in monthly gridded land and marine temperatures
;
print,'Reading gridded temperature'
ncid=ncdf_open(fn1)
fd40=crunc_rddata(ncid,tst=[yrstart,0],tend=[yrend,11],grid=g)  
ncdf_close,ncid
;
; Compute JJA average at each box
;
print,'Computing JJA means'
fd40=reform(fd40,g.nx*g.ny,nmon,nyr)
if dojja eq 0 then fdjja=mkseason(fd40,3,8,datathresh=4) $
              else fdjja=mkseason(fd40,5,7,datathresh=2)
fdjja=reform(fdjja,g.nx,g.ny,nyr)
;
; Extract Carribean grid box time series
;
keepi=where(g.x eq -62.5)
keepj=where(g.y eq 17.5)
carts=reform(fdjja(keepi,keepj,*))
;
; Compute Northern Hemisphere mean
;
print,'Computing NH mean temperature'
nhjkeep=where(g.y ge 0.)
nhy=g.y(nhjkeep)
fd1=fdjja(*,nhjkeep,*)           ; extract northern hemisphere
fd2=globalmean(fd1,nhy)         ; average NH, weighted by latitude
nhlmyr=timey
nhlmts=fd2
;
; Instead of using re-scaled NH10, using regression-based version
;
print,correlate(nh10ts,nhlmts)
acoeff=linfit(nh10ts,nhlmts)
nh10ts=acoeff(0)+acoeff(1)*nh10ts
;
; Compute patterns of correlation coefficients
;
fdjja=reform(fdjja,g.nx*g.ny,nyr)
allylat=fltarr(g.nx,g.ny)
for iy = 0 , g.ny-1 do allylat(*,iy)=g.y(iy)
allylat=reform(allylat,g.nx*g.ny)
print,'Computing full correlation patterns'
patterndata,fdjja,nhlmts,zcorr=nhcorr,ylat=allylat
nhcorr=reform(nhcorr,g.nx,g.ny)
print,'Computing residual correlation patterns'
nhrests=nhlmts-nh10ts
patterndata,fdjja,nhrests,zcorr=rescorr,ylat=allylat
rescorr=reform(rescorr,g.nx,g.ny)
;
loadct,39
def_1color,20,color='vlgreen'
def_1color,21,color='red'
def_1color,22,color='deepblue'
multi_plot,nrow=4,layout='large'
if !d.name eq 'X' then begin
  window, ysize=800
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=14
endelse
;
labels=def_labels(/off)
map=def_map() & map.limit=[-90,-170,90,190.] & map.origx=10.
sm=def_sm() & sm.thresh=0.1 & sm.wid=5
coast=def_coast(/get_device) & coast.double=0 & coast.fill=1
coast.fillcolor=20
levs=[-9,-8,-7,-6,-5,-4,-3,-2,-1,1,2,3,4,5,6,7,8,9]*0.1
levt=[replicate(4,8),2,2,replicate(4,8)]
levl=intarr(18)+1
levc=[replicate(22,9),replicate(21,9)]
;
!p.multi(0)=0
plot,[0,1],/nodata,xstyle=4,ystyle=4
plot,nh10yr,nh10ts,$
  /xstyle,xtitle='Year',$
  /ystyle,yrange=[-0.7,0.4],ytitle='Temperature anomaly  (!Uo!NC)',$
  title='Northern Hemisphere '+titadd+' temperature'
oplot,thick=5,nhlmyr,nhlmts
oplot,!x.crange,[0.,0.]
;
!p.multi=[1,1,2,0,0]
inter_confd,nhcorr,g.x,g.y,$
  sm=sm,map=map,labels=labels,coast=coast,$
  levels=levs,/follow,/hi_on,c_thick=levt,c_label=levl,miss_grey='white',$
  c_colors=levc,$
  title='Correlation with '+titadd+' gridded temperatures'
;
if !d.name eq 'PS' then font_size=10
print,'lon'
lon_cylindrical,map=map,findgen(6)*60.-120.,/slabel
print,'lat'
lat_cylindrical,map=map,findgen(7)*30.-90.,/wlabel
if !d.name eq 'PS' then font_size=14
;
!p.multi=[0,1,4,0,0]
pause
plot,[0,1],/nodata,xstyle=4,ystyle=4
plot,nh10yr,nhrests,thick=5,$
  /xstyle,xtitle='Year',$
  /ystyle,yrange=[-0.7,0.4],ytitle='Residual temperature anomaly  (!Uo!NC)',$
  title='Residual Northern Hemisphere '+titadd+' temperature'
oplot,!x.crange,[0.,0.]
; Compare with Carribean time series
kl=where(finite(carts+nhrests),nkeep)
print,correlate(carts(kl),nhrests(kl))
aaa=linfit(carts(kl),nhrests(kl))
;oplot,nhlmyr,aaa(0)+aaa(1)*carts
;
!p.multi=[1,1,2,0,0]
inter_confd,rescorr,g.x,g.y,$
  sm=sm,map=map,labels=labels,coast=coast,$
  levels=levs,/follow,/hi_on,c_thick=levt,c_label=levl,miss_grey='white',$
  c_colors=levc,$
  title='Correlation with '+titadd+' gridded temperatures'
;
if !d.name eq 'PS' then font_size=10
print,'lon'
lon_cylindrical,map=map,findgen(6)*60.-120.,/slabel
print,'lat'
lat_cylindrical,map=map,findgen(7)*30.-90.,/wlabel
if !d.name eq 'PS' then font_size=14
;
end