;
; Plots a HovMueller diagram (longitude-time) of meridionally averaged
; growing season reconstructions. Uses "corrected" MXD - but shouldn't usually
; plot past 1960 because these will be artificially adjusted to look closer to
; the real temperatures.
;
dosmooth=1         ; smooth each longitude in time?
thalf=10
;
doinfill=1         ; use PCR-infilled data or not?
;
donewref=0         ; use a new reference period or not?
newref=[1801,1900]
;
wlim=-180          ; window to use
elim=180
slim=45
nlim=90
;
doinstr=1          ; use real observations or not?
domask=1           ; mask real observations with reconstruction grid?
;
; Now prepare for plotting
;
loadct,39
multi_plot,nrow=1,layout='large'
if !d.name eq 'X' then begin
  window,ysize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=10
endelse
def_1color,20,color='palepurple'
def_1color,21,color='lpurple'
def_1color,22,color='deepblue'
def_1color,23,color='mlblue'
def_1color,24,color='vlblue'
def_1color,25,color='vvlgreen'
def_1color,26,color='lsand'
def_1color,27,color='orange'
def_1color,28,color='red'
def_1color,29,color='dred'
def_1color,40,color='mlgrey'
;
levs=[-100,-2,-1.2,-0.8,-0.4,-0.2,0,0.3,0.6,1,100]
if donewref ne 0 then levs=levs+0.2
cols=indgen(10)+20
cw=40
cb=!p.color
coll=[cw,cw,cw,cw,cb,cb,cb,cb,cb,cw,cw]
;
; Get the calibrated data
;
print,'Reading reconstructions'
if doinfill eq 0 then begin
  restore,'calibmxd5.idlsave'
  ; Gets: g,mxdyear,mxdnyr,fdcalibu,fdcalibc,mxdfd2,timey,fdseas
endif else begin
  restore,'calibmxd5_pcr.idlsave'
  ; Gets: g,mxdyear,mxdnyr,fdcalibc,timey,fdseas
endelse
;
; Use (and optionally mask) the instrumental data if required
;
if doinstr ne 0 then begin
  if domask eq 0 then begin
    fdcalibc=fdseas
    mxdyear=timey
  endif else begin
    kmxd=where(mxdyear ge timey(0),nyrmxd)
    ktem=where(timey le mxdyear(mxdnyr-1),nyrtem)
    if nyrmxd ne nyrtem then message,'Oooops!'
    fdmask=fdcalibc(*,*,kmxd)
    fdcalibc=fdseas(*,*,ktem)
    mxdyear=timey(ktem)
    ml=where(finite(fdmask) eq 0,nmiss)
    if nmiss gt 0 then fdcalibc(ml)=!values.f_nan
  endelse
endif
;
; Now extract the required window
;
print,'Extracting longitude window'
kx=where((g.x ge wlim) and (g.x le elim),nx)
xlon=g.x(kx)
fd=fdcalibc(kx,*,*)
;
print,'Extracting latitude window'
ky=where((g.y ge slim) and (g.y le nlim),ny)
fd=fd(*,ky,*)
print,'Meridional averaging'
fdtot=total(fd,2,/nan)
fdnum=total(finite(fd),2)
fd=fdtot/float(fdnum)
;
; Optionally re-reference the diagram to a new base period
;
if donewref ne 0 then begin
  print,'Re-referencing to a new base period',newref
  for ix = 0 , nx-1 do begin
    print,ix,format='($,I4)'
    xxx=reform(fd(ix,*))
    mkanomaly,xxx,mxdyear,refperiod=newref
    fd(ix,*)=xxx(*)
  endfor
  print
endif
;
; Optionally smooth the HovMueller diagram in the time direction
;
if dosmooth ne 0 then begin
  print,'Smoothing in the time direction'
  for ix = 0 , nx-1 do begin
    print,ix,format='($,I4)'
    xxx=reform(fd(ix,*))
    filter_cru,thalf,/nan,tsin=xxx,tslow=xlo
    fd(ix,*)=xlo(*)
  endfor
  print
endif
;
contour,fd,xlon,mxdyear,$
  /xstyle,xtitle='Longitude',$
  yrange=[1400,1960],/ystyle,ytitle='Year',$
  levels=levs,c_colors=cols,/cell_fill,$
  ymargin=[10,2]
if dosmooth ne 0 then begin
  contour,fd,xlon,mxdyear,/overplot,levels=levs,c_colors=coll
endif
;
if !d.name eq 'PS' then device,font_size=8
pertit='1961-90'
if donewref ne 0 then pertit=string(newref(0),format='(I4)')+'-'+$
                             string(newref(1),format='(I4)')
!p.multi(0)=1
!p.position=[0.3,0.015,0.7,0.06-0.015]
scale_horiz,levels=levs,c_colors=cols,noextremes=['Below','Above'],$
  title='!Uo!NC wrt '+pertit
!p.position=[0,0,0,0]
;
end