;
; 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
dointerp=1         ; optionally interpolate to give finer longitudinal res
;
doinfill=1         ; use PCR-infilled data or not?
doabd=1          ; use ABD-adjusted data or not? 0=no, 1=yes, -1=show ABD minus Hug
docorr=1           ; use corrected version or not? (uncorrected only available
                   ;                                for doinfill=doabd=0)
;
donewref=1         ; use a new reference period or not?
newref=[1801,1900]
;
wlim=-180          ; window to use
elim=180
slim=60
nlim=90
dohalf=0           ; 0=full A4 1=half A4
dolabel=1          ; 0=do not, 1=do label North America, Europe and Russia
;
doinstr=0          ; use real observations or not?
domask=0           ; mask real observations with reconstruction grid? or reconstruction with obs grid
yr=[1855,1995] ; usually use 1400,1960
yr=[1855,1960] ; usually use 1400,1960
yr=[1400,1960] ; usually use 1400,1960
;yr=[1856,1960] ; usually use 1400,1960
;yr=[1856,1910] ; usually use 1400,1960
;yr=[1800,1899] ; usually use 1400,1960
;
docol=1            ; 0=B&W, 1=colour
;
; Now prepare for plotting
;
loadct,39
multi_plot,nrow=1,layout='large'
if !d.name eq 'X' then begin
  wintim,doinstr*2,ysize=750-350*dohalf
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=10
  if dohalf then device,ysize=14
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'
def_1color,41,color='mdgrey'
def_1color,42,color='lgrey'
;
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
if docol eq 0 then begin
  cw=!p.background
  cb=!p.color
  cs=41
  cg=42
  cols=[cw,cw,cw,cw,cw,cs,cs,cs,cs,cs]
  coll=[cb,cb,cb,cb,cb,cb,cg,cg,cg,cg,cg]
endif else begin
  cols=indgen(10)+20
  cw=40
  cb=!p.color
  coll=[cw,cw,cw,cw,cb,cb,cb,cb,cb,cw,cw]
endelse
;
; Get the calibrated data
;
print,'Reading reconstructions'
if doabd le 0 then begin
  if doinfill eq 0 then begin
    restore,'calibmxd5.idlsave'
    ; Gets: g,mxdyear,mxdnyr,fdcalibu,fdcalibc,mxdfd2,timey,fdseas
    if docorr eq 0 then fdcalibc=fdcalibu
  endif else begin
    restore,'calibmxd5_pcr.idlsave'
    ; Gets: g,mxdyear,mxdnyr,fdcalibc,timey,fdseas
  endelse
endif
;
if doabd lt 0 then begin
  hugfd=fdcalibc
  hugyr=mxdyear
  hugnyr=mxdnyr
endif
;
if doabd ne 0 then begin
  if doinfill eq 0 then begin
    restore,'../mann/calibmxd5_abdlow.idlsave'
    ; Gets: g,mxdyear,mxdnyr,fdcalibu
  endif else begin
    restore,'../mann/calibmxd5_abdlow_pcr.idlsave'
    ; Gets: g,mxdyear,mxdnyr,fdcalibu
  endelse
endif
;
if doabd lt 0 then begin
  if hugnyr ne mxdnyr then message,'mismatch1'
  if total(abs(hugyr-mxdyear)) gt 0 then message,'mismatch2'
  fdcalibc=fdcalibc-hugfd
endif
;
; Optionally mask the reconstructed data if required
;
if doinstr eq 0 then begin
  if domask ne 0 then begin
    kmxd=where(mxdyear ge timey(0),nyrmxd)
    ktem=where(timey le mxdyear(mxdnyr-1),nyrtem)
    if nyrmxd ne nyrtem then message,'Oooops!'
    fdmask=fdseas(*,*,ktem)
    fdcalibc=fdcalibc(*,*,kmxd)
    mxdyear=mxdyear[kmxd]
    ml=where(finite(fdmask) eq 0,nmiss)
    if nmiss gt 0 then fdcalibc(ml)=!values.f_nan
    mxdnyr=nyrmxd
  endif
endif
;
; 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
  mxdnyr=n_elements(mxdyear)
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
;
if dointerp ne 0 then begin
  print,'Interpolating to a finer longitudinal resolution'
  ; We just use simple linear interpolation in a longitudinal direction.  The
  ; main reason is so that if we have one time series surrounded by missing
  ; data longitudinally, we can still plot some contours.
  ;
  nxnew=nx*2+1
  iele=indgen(nx)*2+1
  xnew=findgen(nxnew)
  xnew(0)=xlon(0)-2.5
  xnew(iele)=xlon(*)
  xnew(iele+1)=xlon(*)+2.5
  print,xnew,format='(10F8.3)'
  ;
  fdnew=fltarr(nxnew,mxdnyr)
  fdnew(iele,*)=fd(*,*)
  fdnew(0,*)=fd(0,*)
  fdnew(nxnew-1,*)=fd(nx-1,*)
  for ix = 0 , nx-2 do begin
    fd2col=fd(ix:ix+1,*)
    totv=total(fd2col,/nan,1)
    numv=total(finite(fd2col),1)
    totv=totv/float(numv)
    ml=where(numv eq 0,nmiss)
    if nmiss gt 0 then totv(ml)=!values.f_nan
    fdnew(iele(ix)+1,*)=totv(*)
  endfor
  ;
  nx=nxnew
  xlon=xnew
  fd=fdnew
endif
;
if n_elements(keepfd) eq n_elements(fd) then begin
  kper=where((mxdyear ge yr[0]) and (mxdyear le yr[1]),nkper)
  cfd1=keepfd[*,kper]
  cfd2=fd[*,kper]
  kcorr=where(finite(cfd1+cfd2))
  print,correlate(cfd1[kcorr],cfd2[kcorr])
endif
keepfd=fd 
;
contour,fd,xlon,mxdyear,$
  xrange=[wlim,elim],/xstyle,xtitle='Longitude',$
  yrange=yr,/ystyle,ytitle='Year',$
  levels=levs,c_colors=cols,/cell_fill,$
  ymargin=[10,2]
if (dosmooth ne 0) and (thalf gt 5) then begin
;  contour,fd,xlon,mxdyear,/overplot,levels=levs,c_colors=coll
endif
if dolabel ne 0 then begin
  axis,xaxis=1,xticks=2,xtickv=[-100,20,100],$
    xtickname=['North America','Europe','Russia']
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