;
; DOESN'T ACTUALLY SAVE ANY RESULTS, JUST MAKES THE PLOTS!!!!
;
; Reads in the gridded Hugershoff MXD data, plus the regional age-banded and
; regional Hugershoff series and attempts to adapt the gridded Hugershoff
; data to have the same low-frequency variability as the ABD regions.
; The procedure is as follows:
;
; HUGREG=Hugershoff regions, ABDREG=age-banded regions, HUGGRID=Hugershoff grid
; The calibrated (uncorrected) versions of all these data sets are used.
; However, the same adjustment is then applied to the corrected version of
; the grid Hugershoff data, so that both uncorrected and corrected versions
; are available with the appropriate low frequency variability.  There is some
; ambiguity during the modern period here, however, because the corrected
; version has already been artificially adjusted to reproduce the largest
; scales of observed temperature over recent decades - so a new adjustment
; would be unwelcome.  Therefore, the adjustment term is scaled back towards
; zero when being applied to the corrected data set, so that it is linearly
; interpolated from its 1950 value to zero at 1970 and kept at zero thereafter.
;
; (1) Compute regional means of HUGGRID to (hopefully) confirm that they
; give a reasonably good match to HUGREG.  If so, then for the remainder of
; this routine, HUGREG is replaced by the regional means of HUGGRID.
;
; (2) For each region, low-pass filter (30-yr) both HUGREG and ABDREG,
; and difference them.  This is the additional low frequency information
; that the Hugershoff data set is missing.
;
; (3) To each grid box in HUGGRID, add on a Gaussian-distance-weighted
; mean of nearby regional low frequency, assuming that the low frequency
; information obtained from (2) applies to a point central to each region.
;
; (4) Compute regional means of the adjusted HUGGRID and confirm that they
; give a reasonable match to ABDREG.
;
; For some regions (CAS, TIBP) the low frequency signal is set to zero because
; the gridded data gives a quite different time series than either of the
; regional-mean series.  Also, for those series limited by the availability
; of age-banded results, I set all values from 1400 to 50 years prior to the
; first non-missing value to zero, and then linearly interpolate this 50 years
; and any other gaps with missing values.  Any missing values at the end of 
; the series are filled in by repeating the final non-missing value.
;
thalf=30.
dolfplot=0             ; if set to 1 then also plot the low frequency
                       ; adjustments curves that are applied
doadj=0   ; 0=don't 1=do adjust variance for sample when computing area-means
;
doreg=[0,2,7]
ndoreg=n_elements(doreg)
;
pantit='('+['a','b','c','d','e','f']+') '
;
; 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=12
endelse
def_1color,20,color='mred'
def_1color,21,color='mdgreen'
def_1color,22,color='lpurple'
def_1color,23,color='brown'
;
; 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 with ABD added
;
print,'Reading gridded MXD'
restore,filename='calibmxd5_abdlow.idlsave'
;  g,mxdyear,mxdnyr,fdcalibu,fdcalibc,mxdfd2,timey,fdseas
fdabd=fdcalibu
;
; Get the gridded calibrated Hugershoff MXD
;
print,'Reading gridded MXD'
restore,filename='../summer_modes/calibmxd5.idlsave'
;  g,mxdyear,mxdnyr,fdcalibu,fdcalibc,mxdfd2,timey,fdseas
;
; Next get the calibrated regional Hugershoff MXD
;
print,'Reading regional MXD'
restore,filename='../regtemp_calibrated.allversion.idlsave'
; 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)
;
; Get rid of last 6 years (1995-2000), so it stops in 1994
;
hugnyr=nyr-6
hugtimey=timey[0:hugnyr-1]
hugregts=hugregts[0:hugnyr-1,*]
;
; Next get the calibrated regional age-banded MXD
;
print,'Reading age-banded MXD'
restore,filename='../bandtemp_calibrated.idlsave'
; Gets: nyr,nreg,calregts,regname,timey,tempregts,tempnyr,temptimey
;
; Drop 1995 as the other data sets don't have it, and just keep 9 regions
;
nreg=9
regname=regname(0:8)
nyr=nyr-1
timey=timey(0:nyr-1)
calregts=calregts(0:nyr-1,0:8)
;
if (mxdnyr ne hugnyr) or (mxdnyr ne nyr) then message,'Times do not match'
;
; Create a mask of those boxes that have at least some data in them
; (the mask is the same for corrected data too).
; (and also for the ABD-adjusted data too, I think, given that this is all
; pre-PCR-extension)
;
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 also locate the
; centroid of each region.  Also plot and correlate them with the other
; regional series.
;
print,'Computing regional means of gridded data'
gridregts=fltarr(nyr,nreg)
gridregts(*,*)=!values.f_nan
abdgridregts=fltarr(nyr,nreg)
abdgridregts(*,*)=!values.f_nan
allx=fltarr(g.nx,g.ny)
for iy = 0 , g.ny-1 do allx(*,iy)=g.x(*)
ally=fltarr(g.nx,g.ny)
for ix = 0 , g.nx-1 do ally(ix,*)=g.y(*)
xcentroid=fltarr(nreg)
ycentroid=fltarr(nreg)
lowfreqts=fltarr(nyr,nreg)
lowfreqtsc=fltarr(nyr,nreg)       ; corrected version goes to zero in 1970
for jreg = 0 , ndoreg-1 do begin
  ireg=doreg[jreg]
  rname=hugregname(ireg)
  print,regname(ireg),' = ',rname
  ;
  ; Compute regional-mean timeseries
  maskfd=boxregs(*,*,ireg)        ; in the region? 1=yes NaN=no
  if doadj eq 0 then begin
    ; Straightforward area average with latitude weighting
    gmts=globalmean(fdcalibu,g.y,mask=maskfd)
    maskfd=boxregs(*,*,ireg)        ; in the region? 1=yes NaN=no
    abdgmts=globalmean(fdabd,g.y,mask=maskfd)
  endif else begin
    ; Variance-adjusted average, with latitude weighting too
    nummxd=total(finite(fdcalibu),3)
    kbox=where((nummxd gt 0) and finite(maskfd),nbox)
    regy=ally(kbox)
    regbox=fltarr(nbox,nyr)
    for iyr = 0 , nyr-1 do begin
      onefd=reform(fdcalibu(*,*,iyr))
      regbox(*,iyr)=onefd(kbox)
    endfor
    gmts=var_adjust(regbox,timey,weight=cos(regy*!dtor),refperiod=[1881,1960],$
      /no_anomaly,/mkrbar,rbar=rbarvalue)
    gmts=gmts*sqrt(rbarvalue)      ; convert to variance of area-mean
  endelse
  gridregts(*,ireg)=gmts(*)
  abdgridregts(*,ireg)=abdgmts(*)
  ;
  ; When locating centroids, take care that x does not cross dateline in ESIB
  kl=where(finite(maskfd*maskmxd),ngot)
  ceny=total(ally(kl))/float(ngot)
  listx=allx(kl)
  if rname eq 'ESIB' then begin
    lowlist=where(listx lt 0,nlow)
    if nlow gt 0 then listx(lowlist)=listx(lowlist)+360.
  endif
  cenx=total(listx)/float(ngot)
  print,'Centroid of',ngot,' boxes is x y',cenx,ceny
  xcentroid(ireg)=cenx
  ycentroid(ireg)=ceny
  ;
  ; Now plot them, after smoothing
  ;
  x1=gridregts(*,ireg)
  x2=hugregts(*,ireg)
  x3=calregts(*,ireg)
  x4=abdgridregts(*,ireg)
  filter_cru,thalf,/nan,tsin=x1,tslow=xl1
  filter_cru,thalf,/nan,tsin=x2,tslow=xl2
  filter_cru,thalf,/nan,tsin=x3,tslow=xl3
  filter_cru,thalf,/nan,tsin=x4,tslow=xl4
  ;
  xdiff=xl3-xl1
  if (rname eq 'CAS') or (rname eq 'TIBP') then xdiff(*)=0.
  ; Check if we have complete data or not
  kl=where(finite(xdiff),nkeep)
  if nkeep lt nyr-1 then begin
    ; If not, then check if there is a gap at the beginning
    ist=kl(0)
    if ist gt 0 then begin    ; Fill early missing values with zero
      i50=(ist-50) > 0
      xdiff(0:i50)=0.
    endif
    ; Check if there is a gap at the end
    ien=kl(nkeep-1)
    if ien lt nyr-1 then begin     ; Fill final missing values with last value
      xdiff(ien+1:nyr-1)=xdiff(ien)
    endif
    ; Interpolate to fill all missing values
    kl=where(finite(xdiff))
    xold=xdiff(kl)
    yrold=timey(kl)
    xnew=interpol(xold,yrold,timey)
    xdiff=xnew
  endif
  lowfreqts(*,ireg)=xdiff(*)
  ;
  ; Corrected version of adjustment series should go linearly from 1950 value
  ; to zero in 1970 and stay zero thereafter
  ;
  xold=xdiff
  kzero=where(timey ge 1970)
  xold(kzero)=0.
  kl=where((timey le 1950) or (timey ge 1970))
  xnew=interpol(xold(kl),timey(kl),timey)
  lowfreqtsc(*,ireg)=xnew(*)
  ;
  yr=[min([xl1,xl2,xl3,xl4],/nan),max([xl1,xl2,xl3,xl4],/nan)]
  pause
  plot,timey,x1,/nodata,$
    /xstyle,xtitle='Year',$
    /ystyle,yrange=yr,$
    ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
    title=pantit[jreg]+hugregname(ireg)
  oplot,timey,xl1,color=20,thick=3
  oplot,timey,xl2,color=21,thick=3
  oplot,timey,xl3,color=22,thick=3
  oplot,timey,xl4,color=23,thick=3
  oplot,!x.crange,[0,0],linestyle=1
  r12=mkcorrelation(x1,x2)
  r13=mkcorrelation(x1,x3)
  r23=mkcorrelation(x2,x3)
  l12=mkcorrelation(xl1,xl2)
  l13=mkcorrelation(xl1,xl3)
  l23=mkcorrelation(xl2,xl3)
  xstr=string([r12,r13,r23,l12,l13,l23],format='(6F8.2)')
;  xyouts,1450,0.3,'r: '+xstr,charsize=0.5
  legpos=convert_coord([1700],[!y.crange[0]+0.4],/data,/to_normal)
  legend,position=legpos,/horiz,['GridHug','Hug','ABD','GridABD'],thick=[3,3,3,3],color=[20,21,22,23]
  ;
  if 0 eq 1 then begin
  plot,timey,xdiff,thick=5,$
    /xstyle,xtitle='Year',$
    /ystyle,yrange=[-2,1],$
    ytitle='Low-frequency temperature signal (!Uo!NC)',$
    title=hugregname(ireg)
  oplot,timey,xnew,thick=2
  oplot,!x.crange,[0,0],linestyle=1
  endif
  ;
endfor
;
end