;
; Computes the annual range (Apr-Sep mean minus Oct-Mar mean) and plots
; timeseries and maps as necessary
;
doinfill=0    ; 0=no infilling allowed, 1=use data with some spatio-temp infil
domap=1       ; 0=don't plot the maps 1=do
donorm=0      ; 0=don't 1=do normalise over 1961-90 prior to differencing
domask=0   ; 0=do no mask, >0=mask with boxes with >=5yr data in [X,X+9]
doland=1     ; -1=ocean 0=all 1=land
don20=1      ; 0=>0N, 1=>20N
;
doold=2        ; 0=use new data set (2002), 1=use old version, 2=use land-only
case doold of
  0: fnadd=''
  1: fnadd='..old..'
  2: fnadd='..land..'
endcase
;
restore,filename='obs_temp_as.idlsave'+fnadd
; Gets: timey,fdseas,xlon,ylat,nyr,nx,ny,$
;  fdltm,fdsd,landonly,seaonly,mainlysea,missfrac,rawseas
if doinfill eq 0 then fdseas=rawseas
;
fdas=fdseas
if donorm eq 1 then begin
  fdas=reform(fdas,nx*ny,nyr)
  mknormal,fdas,timey,refperiod=[1961,1990]
  fdas=reform(fdas,nx,ny,nyr)
endif
;
restore,filename='obs_temp_om.idlsave'+fnadd
if doinfill eq 0 then fdseas=rawseas
;
fdom=fdseas
if donorm eq 1 then begin
  fdom=reform(fdom,nx*ny,nyr)
  mknormal,fdom,timey,refperiod=[1961,1990]
  fdom=reform(fdom,nx,ny,nyr)
endif
;
fdar=fdas-fdom
;
; Find mask if required
;
if domask gt 0 then begin
  masklist=where((timey ge domask) and (timey le domask+9),nmask)
  if nmask ne 10 then message,'Ooops!'
  fdmask=fdar(*,*,masklist)
  masknum=total(finite(fdmask),3)
  ml=where(masknum lt 5,nmiss)
  fdmask=fltarr(nx,ny)
  if nmiss gt 0 then fdmask(ml)=!values.f_nan
endif
;
; Compute hemispheric means
;
if don20 eq 0 then ky=where(ylat ge 0.) else ky=where(ylat ge 20.)
;
if (domask gt 0) or (doland ne 0) then begin
  if (domask eq 0) then begin
    if doland lt 0 then fdmask=seaonly else fdmask=landonly
  endif else begin
    if doland lt 0 then fdmask=fdmask*seaonly
    if doland gt 0 then fdmask=fdmask*landonly
  endelse
  fd=fdmask
  dummy=globalmean(fdas[*,ky,*],ylat[ky],nhemi=nhas,mask=fd[*,ky])
  fd=fdmask
  dummy=globalmean(fdom[*,ky,*],ylat[ky],nhemi=nhom,mask=fd[*,ky])
  fd=fdmask
  dummy=globalmean(fdar[*,ky,*],ylat[ky],nhemi=nhar,mask=fd[*,ky])
endif else begin
  dummy=globalmean(fdas[*,ky,*],ylat[ky],nhemi=nhas)
  dummy=globalmean(fdom[*,ky,*],ylat[ky],nhemi=nhom)
  dummy=globalmean(fdar[*,ky,*],ylat[ky],nhemi=nhar)
endelse
;
; Now prepare for plotting
;
loadct,39
multi_plot,nrow=3,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=12
endelse
def_1color,20,color='vlpurple'
def_1color,22,color='purple'
def_1color,24,color='deepblue'
def_1color,27,color='vlblue'
def_1color,28,color='lgrey'
def_1color,29,color='lyellow'
def_1color,31,color='orange'
def_1color,34,color='dred'
def_1color,36,color='lred'
def_1color,100,color='red'
def_smearcolor,fromto=[20,22]
def_smearcolor,fromto=[22,24]
def_smearcolor,fromto=[24,27]
def_smearcolor,fromto=[29,31]
def_smearcolor,fromto=[31,34]
def_smearcolor,fromto=[34,36]
cols=indgen(17)+20
;
if donorm eq 0 then yr=[-1,0.7] else yr=[-1.6,1.6]
if doland gt 0 then yr=yr*2.
if domask gt 0 then mtit='Masked with coverage from '+string(domask,$
  format='(I4)')+'s' else mtit=' '
if doland lt 0 then mtit='Ocean-only '+mtit
if doland gt 0 then mtit='Land-only '+mtit
if don20 gt 0 then mtit='>20!Uo!NN '+mtit
plot,timey,nhom,$
  /xstyle,xtitle='Year',$
  /ystyle,ytitle='NH Apr-Sep or Oct-Mar temp  (!Uo!NC wrt 61-90)',$
  yrange=yr,thick=2,title=mtit
oplot,timey,nhas,thick=2,color=100
oplot,!x.crange,[0,0],linestyle=1
filter_cru,10,/nan,tsin=nhom,tslow=tsl
oplot,timey,tsl,thick=6
filter_cru,10,/nan,tsin=nhas,tslow=tsl
oplot,timey,tsl,thick=6,color=100
legend,['Oct-Mar','Apr-Sep'],thick=[2,2],color=[!p.color,100]
;
if donorm eq 0 then yr=[-0.6,0.8] else yr=[-0.8,1.2]
if doland gt 0 then yr=yr*2.
plot,timey,nhar,$
  /xstyle,xtitle='Year',$
  /ystyle,ytitle='NH Apr-Sep minus Oct-Mar temp  (!Uo!NC)',$
  yrange=yr,thick=2,title=mtit
oplot,!x.crange,[0,0],linestyle=1
filter_cru,10,/nan,tsin=nhar,tslow=tsl
oplot,timey,tsl,thick=6
;
; Compute and then plot decadal mean maps of the temperature difference
;
if domap ne 0 then begin
kl=where((timey ge 1860) and (timey le 1999),nyr)
fdar=fdar(*,*,kl)
ndec=nyr/10
fdar=reform(fdar,nx,ny,10,ndec)
fdtot=total(fdar,/nan,3)
fdnum=total(finite(fdar),3)
ml=where(fdnum lt 5,nmiss)
fddec=fdtot/float(fdnum)
if nmiss gt 0 then fddec(ml)=!values.f_nan
;
!p.multi=[0,2,4,0,0]
;levs=findgen(21)*0.5-3.8
levs=[-7,-4,-2,-1.5,-1,-0.8,-0.6,-0.4,-0.2]
levs=[levs,-reverse(levs)]
;if donorm then levs=levs-1.2
map=def_map()
map.limit=[0.,-180,90,180]
map.origx=0.
map.xmargin=[1,1]  & map.ymargin=[1,2]
coast=def_coast(/get_device)  & coast.double=0
labels=def_labels(/off)
;
for idec = 0 , ndec-1 do begin
  pause
  inter_boxfd,fddec(*,*,idec),xlon,ylat,$
    map=map,coast=coast,labels=labels,$
    /scale,levels=levs,c_colors=cols,$
    title=string(idec*10+1860,format='(I4)')+'s'
endfor
endif
;
end