;
; 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
;
; 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'
;
; 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'
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).
;
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
gridregts=fltarr(nyr,nreg)
gridregts(*,*)=!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 ireg = 0 , nreg-1 do begin
  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)
  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(*)
  ;
  ; 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)
  filter_cru,thalf,/nan,tsin=x1,tslow=xl1
  filter_cru,thalf,/nan,tsin=x2,tslow=xl2
  filter_cru,thalf,/nan,tsin=x3,tslow=xl3
  ;
  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(*)
  ;
  pause
  plot,timey,x1,/nodata,$
    /xstyle,xtitle='Year',$
    /ystyle,yrange=[-2,1],$
    ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
    title=hugregname(ireg)
  oplot,timey,xl1,color=20,thick=3
  oplot,timey,xl2,color=21,thick=3
  oplot,timey,xl3,color=22,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
  legend,/horiz,['Grid','Hug','ABD'],thick=[3,3,3],color=[20,21,22]
  ;
  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
  ;
endfor
;
; Now add appropriate low frequency average to each grid point.  The width of
; the Gaussian weighting is selected somewhat arbitrarily, though noting that
; the furthest distance between a grid box and its nearest regional centroid
; varies from 0 to 2000 km.
; A slight complication is that we do not want the CAS and TIBP zero
; adjustments to influence other regions (notably NSIB).  We can't just
; set them to missing or remove them, because then CAS and TIBP boxes would
; be adjusted using the nearest regions which would then be NSIB.  So we need
; to check for the special case where the NEAREST region is either CAS or
; TIBP and then use all regions, OTHERWISE remove CAS and TIBP prior to
; doing the calculation.
;
regcas=where(hugregname eq 'CAS') & regcas=regcas(0)
regtibp=where(hugregname eq 'TIBP') & regtibp=regtibp(0)
regother=where((hugregname ne 'CAS') and (hugregname ne 'TIBP'))
gwid=1200
fdabd=fdcalibu
fdabdc=fdcalibc
multi_plot,nrow=8,ncol=2
!y.margin=[2,1]
for ix = 0 , g.nx-1 do begin
  for iy = 0 , g.ny-1 do begin
    if finite(maskmxd(ix,iy)) then begin
      ; Compute distance (km) between box and all region centroids
      listdist=geog_dist(g.x(ix),g.y(iy),xcentroid,ycentroid)
      listwt=exp(-((listdist/gwid)^2)/2.)/sqrt(2*!pi)
      sumwt=total(listwt)
      listwt=listwt/sumwt
      ; Check whether CAS or TIBP are the closest; if so use all regions
      ; otherwise remove CAS and TIBP from the analysis.
      dummy=min(listdist,minele)
      if (minele eq regcas) or (minele eq regtibp) then begin
        lfts=lowfreqts               ; use all time series
        lftsc=lowfreqtsc             ; use all time series (corrected version)
        regwt=listwt                 ; use all weights
        print,'FULL ',format='($,A5)'
      endif else begin
        lfts=lowfreqts(*,regother)   ; use all years, but only 'other' regions
        lftsc=lowfreqtsc(*,regother) ; use all years, but only 'other' regions
        regwt=listwt(regother)       ; use weights only from 'other' regions
        print,'SOME ',format='($,A5)'
      endelse
      lfreq=regwt##lfts
      lfreqc=regwt##lftsc
      fdabd(ix,iy,*)=fdcalibu(ix,iy,*)+lfreq(*)
      fdabdc(ix,iy,*)=fdcalibc(ix,iy,*)+lfreqc(*)
      x=reverse(listwt(sort(listwt)))
      print,x/x(0),format='(9F8.3)'
      if dolfplot then begin
        pause
        plot,timey,lfreq,/xstyle,title=string(g.x(ix),g.y(iy)),thick=3
        oplot,timey,lfreqc
        oplot,!x.crange,[0,0],linestyle=1
      endif
    endif
  endfor
endfor
!y.margin=[4,2]
;
; Now compute regional means of the new adjusted gridded MXD data
;
print,'Computing regional means of gridded data'
multi_plot,nrow=4,layout='large'
for ireg = 0 , nreg-1 do begin
  rname=hugregname(ireg)
  print,regname(ireg),' = ',rname
  ;
  ; Compute regional-mean timeseries
  maskfd=boxregs(*,*,ireg)        ; in the region? 1=yes NaN=no
  gmts=globalmean(fdabd,g.y,mask=maskfd)
  maskfd=boxregs(*,*,ireg)        ; in the region? 1=yes NaN=no
  gmtsc=globalmean(fdabdc,g.y,mask=maskfd)
  ;
  ; Now plot them, after smoothing
  ;
  x1=gmts
  x1c=gmtsc
  x2=hugregts(*,ireg)
  x3=calregts(*,ireg)
  filter_cru,thalf,/nan,tsin=x1,tslow=xl1
  filter_cru,thalf,/nan,tsin=x1c,tslow=xl1c
  filter_cru,thalf,/nan,tsin=x2,tslow=xl2
  filter_cru,thalf,/nan,tsin=x3,tslow=xl3
  ;
  pause
  plot,timey,x1,/nodata,$
    /xstyle,xtitle='Year',$
    /ystyle,yrange=[-2,1],$
    ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
    title=hugregname(ireg)
  oplot,timey,xl1,color=20,thick=3
  oplot,timey,xl1c,color=20,thick=1
  oplot,timey,xl2,color=21,thick=3
  oplot,timey,xl3,color=22,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
  legend,/horiz,['GridABD','Hug','ABD'],thick=[3,3,3],color=[20,21,22]
  ;
endfor
;
fdcalibu=fdabd
fdcalibc=fdabdc
save,filename='calibmxd5_abdlow.idlsave',$
  g,mxdyear,mxdnyr,fdcalibu,fdcalibc
;
end