; We have previously (calibrate_mxd.pro) calibrated the high-pass filtered
; MXD over 1911-1990, applied the calibration to unfiltered MXD data (which
; gives a zero mean over 1881-1960) after extending the calibration to boxes
; without temperature data (pl_calibmxd1.pro).  We have identified and
; artificially removed (i.e. corrected) the decline in this calibrated
; data set.  We now recalibrate this corrected calibrated dataset against
; the unfiltered 1911-1990 temperature data, and apply the same calibration
; to the corrected and uncorrected calibrated MXD data.
;
calper=[1911,1990]
verper=[1856,1910]
;
; Get corrected and uncorrected calibrated MXD data
;
restore,filename='calibmxd3.idlsave'
;  g,mxdyear,mxdnyr,fdcalib,mxdfd2,fdcorrect
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=fltarr(g.nx,g.ny) & fdcltm(*,*)=!values.f_nan
fdcalibu=fltarr(g.nx,g.ny,mxdnyr) & fdcalibu(*,*,*)=!values.f_nan
fdcalibc=fltarr(g.nx,g.ny,mxdnyr) & fdcalibc(*,*,*)=!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(fdcorrect(ix,iy,*))
      onemxdu=reform(fdcalib(ix,iy,*))      ; uncorrected calibrated MXD
      onetem=reform(fdseas(ix,iy,*))
      mxdhi=onemxd
      temhi=onetem
      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
        fdcorr(ix,iy)=fc1(0)
        fdalph(ix,iy)=fc1(1)
        fdbeta(ix,iy)=fc1(2)
        ltm=total(temhi)/float(nkeep)
        fdcltm(ix,iy)=ltm
        mxdpred=fc1(1)+fc1(2)*mxdhi
        sqvfull=total((temhi-ltm)^2)
        sqvresid=total((temhi-mxdpred)^2)
        fdvexp(ix,iy)=1.-sqvresid/sqvfull
        fdcalibc(ix,iy,*)=fc1(1)+fc1(2)*onemxd(*)
        fdcalibu(ix,iy,*)=fc1(1)+fc1(2)*onemxdu(*)
      endif
    endif
  endfor
endfor
print
;
; Now verify on a grid-box basis
; No need to verify the correct and uncorrected versions, since these
; should be identical prior to 1920 or 1930 or whenever the decline
; was corrected onwards from.
;
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(fdcalibc(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)
        mxdpred=fdalph(ix,iy)+fdbeta(ix,iy)*mxdhi
        sqvfull=total((temhi-ltm)^2)
        sqvresid=total((temhi-mxdpred)^2)
        fdvver(ix,iy)=1.-sqvresid/sqvfull
      endif
    endif
  endfor
endfor
print
;
; Now save the data for later analysis
;
save,filename='calibmxd4.idlsave',$
  g,mxdyear,mxdnyr,fdcorr,fdalph,fdbeta,fdvexp,fdcalibc,fdcalibu,mxdfd2,$
  fdrver,fdvver,fdseas,timey
;
end