;
; Reads in the gridded Hugershoff MXD data, plus the
; regional Hugershoff series and compares them at the regional scale.
;
thalf=30.
matchvar=1        ; use gridded data calibrated using 0=regression 1=variance
;
doabd=1    ; 0=compare with HUG, 1=with ABD
dodiff=1
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=4,layout='large'
if !d.name eq 'X' then begin
  wintim,ysize=800
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=20
endelse
def_1color,20,color='mred'
def_1color,21,color='mdgreen'
def_1color,22,color='lpurple'
;
; Read in pre-computed masks of those grid boxes in each region
;
restore,'../reg_mxdboxes.idlsave'
; Gets: nreg,boxlists,boxnumbs,regname,g,boxregs,boxprec,boxtemp
knh=where(g.y gt 0)
boxregs=boxregs(*,knh,*)
boxlists=boxlists(*,knh,*)
;
; Get the gridded calibrated Hugershoff MXD
;
print,'Reading gridded MXD'
if matchvar eq 0 then fnadd='_regress' else fnadd=''
restore,filename='../summer_modes/calibmxd5'+fnadd+'.idlsave'
;  g,mxdyear,mxdnyr,fdcalibu,fdcalibc,mxdfd2,timey,fdseas
;
; Next get the calibrated regional Hugershoff MXD
;
print,'Reading regional MXD'
if doabd then begin
  restore,filename='../bandtemp_calibrated.idlsave'
endif else begin
  restore,filename='../regtemp_calibrated.allversion.idlsave'
endelse
; Gets: nyr,nreg,calregts,regname,timey,tempregts,tempnyr,temptimey
;
; Keep just the 9 individual regions
;
hugnreg=9
hugregname=regname(0:8)
hugregts=calregts(*,0:8)
hugnyr=nyr
hugtimey=timey
;
; Drop everything after 1994 as they're not in the gridded data
;
kl=where(hugtimey lt 1995,hugnyr)
hugtimey=hugtimey[kl]
hugregts=hugregts[kl,*]
;
if mxdnyr ne hugnyr then message,'Times do no match'
nyr=hugnyr
timey=hugtimey
nreg=hugnreg
;
; Create a mask of those boxes that have at least some data in them
; (the mask is the same for corrected data too).
;
maskmxd=total(finite(fdcalibu),3)
ml=where(maskmxd eq 0)
maskmxd(*,*)=1
maskmxd(ml)=!values.f_nan
;
; Now compute regional means of the gridded MXD data and
; plot them with the other regional series.
;
print,'Computing regional means of gridded data'
gridregts=fltarr(nyr,nreg)
gridregts(*,*)=!values.f_nan
;nreg=2
;for ireg = 0 , nreg-1 do begin
for ireg = 1 , 1 do begin
  rname=hugregname(ireg)
  print,regname(ireg),' = ',rname
  ;
  ; Compute regional-mean timeseries
  maskfd=boxregs(*,*,ireg)        ; in the region? 1=yes NaN=no
  ; Straightforward area average with latitude weighting
  gmts=globalmean(fdcalibu,g.y,mask=maskfd)
  gridregts(*,ireg)=gmts(*)
  ;
  ; Now plot them, after smoothing
  ;
  x1=gridregts(*,ireg)
  x2=hugregts(*,ireg)
;  mkanomaly,x1,timey,refperiod=[1881,1960]
;  mkanomaly,x2,timey,refperiod=[1881,1960]
  filter_cru,thalf,/nan,tsin=x1,tslow=xl1
  filter_cru,thalf,/nan,tsin=x2,tslow=xl2
  ;
  pause
  plot,timey,x1,/nodata,$
    /xstyle,xtitle='Year',$
    /ystyle,yrange=[-1.5,0.5],$
    ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
    title=hugregname(ireg)
  oplot,timey,xl1,color=20,thick=5
  oplot,timey,xl2,color=21,thick=5
  oplot,!x.crange,[0,0],linestyle=1
  legtxt=['Grid','Hug']
  if doabd then legtxt[1]='ABD'
  legend,/horiz,legtxt,thick=[5,5],color=[20,21]
  ;
  if dodiff then begin
    pause
    plot,timey,xl2-xl1,thick=5,$
      /xstyle,xtitle='Year',$
      /ystyle,yrange=[-1.,1.],$
      ytitle='Difference  (!Uo!NC)'
    oplot,!x.crange,[0,0],linestyle=1
  endif
  ;
endfor
;
end