;
; Calibrates the gridded and infilled MXD data against instrumental
; summer temperatures (land&sea).  On a grid-box basis first, using the
; period 1911-1990 for calibration and the period 1856-1910 for verification,
; where data is available.
;
; Due to the decline, all time series are first high-pass filter with a
; 40-yr filter, although the calibration equation is then applied to raw
; data.
;
matchvar=1        ; 0=regression, 1=variance matching
thalf=40.
dofilt=1          ; 0=don't high-pass filter, 1=high-pass filter
calper=[1911,1990]
verper=[1856,1910]
;
loadct,39
multi_plot,nrow=4,ncol=2,layout='large'
if !d.name ne 'PS' then window,ysize=750
;
; Get the density data
;
print,'Restoring MXD data'
restore,'mxd_infill.idlsave'
; Gets: mxdfd2,mxdyear,g and others
nvalues=total(finite(mxdfd2),3)
;
; Get the temperature data
;
print,'Restoring temperature data'
restore,'obs_temp_as.idlsave'
; Gets: timey,fdseas,xlon,ylat,nyri,fdltm and others
; 
; Calibrate grid-box by grid-box
;
print,'Calibrating each grid box',calper
kcalt=where( (timey ge calper(0)) and (timey le calper(1)) , nct )
kcalm=where( (mxdyear ge calper(0)) and (mxdyear le calper(1)) , ncm )
if nct ne ncm then message,'Ooops!'
;
fdcorr=fltarr(g.nx,g.ny) & fdcorr(*,*)=!values.f_nan
fdalph=fltarr(g.nx,g.ny) & fdalph(*,*)=!values.f_nan
fdbeta=fltarr(g.nx,g.ny) & fdbeta(*,*)=!values.f_nan
fdvexp=fltarr(g.nx,g.ny) & fdvexp(*,*)=!values.f_nan
; fdcltm is the mean temperature over the calibration period, required for
; later calculation of the RE, and also for working out the additive
; factor needed to make the reconstructed and observed means match over
; the calibration period.
fdcltm=fltarr(g.nx,g.ny) & fdcltm(*,*)=!values.f_nan
; fdcaloffset is the additive offset needed to set the reconstruction
; calibration period mean equal to, so that the reconstruction has the
; same mean as the observed temperature during their overlap within the
; calibration period.  fdcaloffset is not the same as fdcltm, because
; missing data within the calibration period means that the calibration is
; not actually done over the full 1911-1990 period.
fdcaloffset=fltarr(g.nx,g.ny) & fdcaloffset(*,*)=!values.f_nan
; fdcalib is the calibration applied to the full, unfiltered MXD
fdcalib=fltarr(g.nx,g.ny,mxdnyr) & fdcalib(*,*,*)=!values.f_nan
;
i=0
for ix = 0 , g.nx-1 do begin
  for iy = 0 , g.ny-1 do begin
    if finite(fdltm(ix,iy)) and (nvalues(ix,iy) gt 20.) then begin
      onemxd=reform(mxdfd2(ix,iy,*))
      onetem=reform(fdseas(ix,iy,*))
      ;
      ; Let's just do a visual check of the temperature series, in case
      ; the early values are weird!
      ;
      pause
      plot,timey,onetem,/xstyle,psym=10,$
        title=string([xlon[ix],ylat[iy]],format='(F7.2," Lon  ",F7.2," Lat")')
      oplot,timey,onetem,psym=1
      oplot,!x.crange,[0,0],linestyle=1
      cutoff=0
      ; Delete obviously erroneous early values at two W US boxes
      if (xlon[ix] eq -117.5) and (ylat[iy] eq 42.5) then cutoff=1874
      if (xlon[ix] eq -107.5) and (ylat[iy] eq 37.5) then cutoff=1871
      if cutoff gt 0 then begin
        iend=where(timey eq cutoff)
        iend=iend[0]
        onetem[0:iend]=!values.f_nan
        oplot,timey,onetem,color=240,psym=10
        oplot,timey,onetem,color=240,psym=1
        xyouts,1900,0,'CHECK THIS. DROPPED ALL BUT THE RED',color=240
        fdseas[ix,iy,0:iend]=!values.f_nan
      endif
      ;
      if dofilt eq 1 then begin
        filter_cru,thalf,/nan,tsin=onemxd,tshigh=mxdhi
        filter_cru,thalf,/nan,tsin=onetem,tshigh=temhi
      endif else begin
        mxdhi=onemxd
        temhi=onetem
      endelse
      mxdraw=onemxd[kcalm]
      temraw=onetem[kcalt]
      mxdhi=mxdhi(kcalm)
      temhi=temhi(kcalt)
      kl=where(finite(mxdhi+temhi),nkeep)
      if nkeep gt 20. then begin
        i=i+1 & print,i,format='($,I4)'
        mxdhi=mxdhi(kl)
        temhi=temhi(kl)
        mkcalibrate,mxdhi,temhi,fitcoeff=fc1,autocoeff=ac1,matchvar=matchvar
        fdcorr(ix,iy)=fc1(0)
        fdalph(ix,iy)=fc1(1)
        fdbeta(ix,iy)=fc1(2)
        ltm=total(temhi)/float(nkeep)
        mxdpred=fc1(1)+fc1(2)*mxdhi
        sqvfull=total((temhi-ltm)^2)
        sqvresid=total((temhi-mxdpred)^2)
        fdvexp(ix,iy)=1.-sqvresid/sqvfull
        ; Store the mean of the unfiltered temperature over the cal period
        fdcltm(ix,iy)=mean(temraw[kl])
        ; Now apply the calibration to the unfiltered series
        fdcalib(ix,iy,*)=fc1(1)+fc1(2)*onemxd(*)
        ; Anomalise the reconstructed against the full calibration period
        ; (excluding missing MXD values, but not excluding missing temperature
        ; values), then compute the mean of the anomalised reconstruction
        ; over the actual calibration subset.  Use this to work out what
        ; offset should be applied to the calibrated values to give the
        ; correct mean level.
        onefullrecon=reform(fdcalib[ix,iy,*])
        mkanomaly,onefullrecon,mxdyear,refperiod=calper,refmean=rfm
        onesubsetrecon=fc1[1]+fc1[2]*mxdraw[kl]
        fdcaloffset[ix,iy]=fdcltm[ix,iy]-mean(onesubsetrecon-rfm[0])
      endif
    endif
  endfor
endfor
print
;
; Now verify on a grid-box basis
;
print,'Verifying each grid box',verper
kcalt=where( (timey ge verper(0)) and (timey le verper(1)) , nct )
kcalm=where( (mxdyear ge verper(0)) and (mxdyear le verper(1)) , ncm )
if nct ne ncm then message,'Ooops!'
fdrver=fltarr(g.nx,g.ny) & fdrver(*,*)=!values.f_nan
fdvver=fltarr(g.nx,g.ny) & fdvver(*,*)=!values.f_nan
i=0
for ix = 0 , g.nx-1 do begin
  for iy = 0 , g.ny-1 do begin
    if finite(fdbeta(ix,iy)) then begin
      onemxd=reform(fdcalib(ix,iy,*))
      ; But the calibrated series had the high-pass calibration applied to
      ; the raw MXD, which would have left the reconstruction having a
      ; near zero mean over 1881-1960, while we would prefer to match the
      ; observed temperature mean over the calibration period.  So let's
      ; impose that.
      mkanomaly,onemxd,mxdyear,refperiod=calper
      onemxd=onemxd+fdcaloffset[ix,iy]
      ;
      onetem=reform(fdseas(ix,iy,*))
      ;
      mxdhi=onemxd(kcalm)
      temhi=onetem(kcalt)
      ;
      kl=where(finite(mxdhi+temhi),nkeep)
      if nkeep gt 9. then begin
        i=i+1 & print,i,format='($,I4)'
        mxdhi=mxdhi(kl)
        temhi=temhi(kl)
        fdrver(ix,iy)=mkcorrelation(mxdhi,temhi)
        ;
        ; Compute RE using calibrated series
        mxdpred=mxdhi
        ltm=fdcltm(ix,iy)    ; must use calibration period mean for RE
        pause
        plot,temhi,/xstyle,$
          title=string([xlon[ix],ylat[iy]],format='(F7.2," Lon  ",F7.2," Lat")')
        oplot,mxdpred,color=240
        ;
        sqvfull=total((temhi-ltm)^2)
        sqvresid=total((temhi-mxdpred)^2)
        fdvver1=1.-sqvresid/sqvfull
        fdvver(ix,iy)=fdvver1
      endif
    endif
  endfor
endfor
print
;
; Now save the data for later analysis
;
if matchvar eq 0 then fnadd='_regress' else fnadd=''
save,filename='calibmxd1'+fnadd+'.idlsave',$
  g,mxdyear,mxdnyr,fdcorr,fdalph,fdbeta,fdvexp,fdcalib,mxdfd2,$
  fdrver,fdvver,fdseas,timey,fdcltm,fdcaloffset,calper
;
end