;
; Computes spatial correlation results (i.e., between gridded MXD series),
; as a function of distance, timescale and local skill.
;
dosmooth=0         ; smooth each box series in time?
thalf=20
;
doinfill=0         ; use PCR-infilled data or not?
doabd=0            ; use ABD-adjusted data or not?
docorr=0           ; use corrected version or not? (uncorrected only available
                   ;                                for doinfill=doabd=0)
;
wlim=-180          ; window to use
elim=180
slim=0
nlim=90
;
doinstr=-1         ; use real observations or not? -ve attempts to factor out T from MXD!
domask=1           ; mask real observations with reconstruction grid?
yr=[1400,2000] ; usually use 1400,1960
yr=[1400,1960] ; usually use 1400,1960
;
; Now prepare for plotting
;
loadct,39
;multi_plot,nrow=1,/landscape
multi_plot,nrow=2
if !d.name eq 'X' then begin
  wintim,xsize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=14
endelse
def_1color,18,color='dgreen'
def_1color,19,color='black'
def_1color,20,color='red'
def_1color,21,color='lpurple'
def_1color,22,color='deepblue'
def_1color,23,color='mlblue'
def_1color,24,color='vlblue'
def_1color,25,color='vvlgreen'
def_1color,26,color='lsand'
def_1color,27,color='orange'
def_1color,28,color='red'
def_1color,29,color='dred'
def_1color,40,color='mlgrey'
;
levs=[-100,-2,-1.2,-0.8,-0.4,-0.2,0,0.3,0.6,1,100]
cols=indgen(10)+20
cw=40
cb=!p.color
coll=[cw,cw,cw,cw,cb,cb,cb,cb,cb,cw,cw]
;
; Get the calibrated data
;
print,'Reading reconstructions'
if doabd eq 0 then begin
  if doinfill eq 0 then begin
    restore,'calibmxd5.idlsave'
    ; Gets: g,mxdyear,mxdnyr,fdcalibu,fdcalibc,mxdfd2,timey,fdseas
    if docorr eq 0 then fdcalibc=fdcalibu
  endif else begin
    restore,'calibmxd5_pcr.idlsave'
    ; Gets: g,mxdyear,mxdnyr,fdcalibc,timey,fdseas
  endelse
endif else begin
  if doinfill eq 0 then begin
    restore,'../mann/calibmxd5_abdlow.idlsave'
    ; Gets: g,mxdyear,mxdnyr,fdcalibu
  endif else begin
    restore,'../mann/calibmxd5_abdlow_pcr.idlsave'
    ; Gets: g,mxdyear,mxdnyr,fdcalibu
  endelse
endelse
;
; Use (and optionally mask) the instrumental data if required
;
if doinstr ne 0 then begin
  print,'Using instrumental data'
  if (domask eq 0) and (doinstr gt 0) then begin
    fdcalibc=fdseas
    mxdyear=timey
  endif else begin
    kmxd=where(mxdyear ge timey(0),nyrmxd)
    ktem=where(timey le mxdyear(mxdnyr-1),nyrtem)
    if nyrmxd ne nyrtem then message,'Oooops!'
    fdmask=fdcalibc(*,*,kmxd)
pause
qtv,total(finite(fdmask),3)
    fdcalibc=fdseas(*,*,ktem)
pause
qtv,total(finite(fdcalibc),3)
    fdcalibc=fdseas(*,*,ktem)
    mxdyear=timey(ktem)
    if doinstr lt 0 then begin
      for ix = 0 , g.nx-1 do begin
        for iy = 0 , g.ny-1 do begin
          ts1=reform(fdmask[ix,iy,*])
          ts2=reform(fdcalibc[ix,iy,*])
          kts=where(finite(ts1+ts2),nkts)
          if nkts ge 20 then begin
            acoeff=linfit(ts2[kts],ts1[kts])
            ts1=ts1-acoeff[0]-acoeff[1]*ts2
            fdcalibc[ix,iy,*]=ts1[*]
          endif else begin
            fdcalibc[ix,iy,*]=!values.f_nan
          endelse
        endfor
      endfor
pause
qtv,total(finite(fdcalibc),3)
    endif else begin
      ml=where(finite(fdmask) eq 0,nmiss)
      if nmiss gt 0 then fdcalibc(ml)=!values.f_nan
    endelse
  endelse
  mxdnyr=n_elements(mxdyear)
endif
;
; Now extract the required window
;
print,'Extracting longitude window'
kx=where((g.x ge wlim) and (g.x le elim),nx)
xlon=g.x(kx)
fd=fdcalibc(kx,*,*)
;
print,'Extracting latitude window'
ky=where((g.y ge slim) and (g.y le nlim),ny)
ylat=g.y[ky]
fd=fd(*,ky,*)
;
; Optionally smooth the series in the time direction
;
if dosmooth ne 0 then begin
  print,'Smoothing in the time direction'
  for ix = 0 , nx-1 do begin
    print,ix,format='($,I4)'
    for iy = 0 , ny-1 do begin
      xxx=reform(fd(ix,iy,*))
      filter_cru,thalf,/nan,tsin=xxx,tslow=xlo
      fd(ix,iy,*)=xlo(*)
    endfor
  endfor
  print
endif
;
; Now compute correlations between available boxes
;
allr=fltarr(100000)
allr[*]=!values.f_nan
alld=fltarr(100000)
nall=0
for ix = 0 , nx-1 do begin
  for iy = 0 , ny-1 do begin
    ts1=fd[ix,iy,*]
    kl=where(finite(ts1),nkeep)
    if nkeep ge 30. then begin
      print,xlon[ix],ylat[iy],nall
      for jx = 0 , nx-1 do begin
        for jy = 0 , ny-1 do begin
          ts2=fd[jx,jy,*]
          kl=where(finite(ts1+ts2),nkeep)
          if nkeep ge 30. then begin
            allr[nall]=correlate(ts1[kl],ts2[kl])
      alld[nall]=geogdist(pos1=[xlon[ix],ylat[iy]],pos2=[xlon[jx],ylat[jy]])
            nall=nall+1
          endif
        endfor
      endfor
    endif
  endfor
endfor
;
; Now fit an overall correlation decay length
;
x=double(alld[0:nall-1])
y=double(allr[0:nall-1])
alow=double(1./15000.)
ahigh=double(1./100.)
iterfit,x,y,alow,ahigh,a
d0=1./a
;
plot,x,y,psym=def_sym(10),symsize=0.3,$
  /xstyle,xrange=[0,5000],xtitle='Separation distance  (km)',$
  /ystyle,yrange=[-0.5,1],ytitle='Correlation'
oplot,!x.crange,[0,0],linestyle=1
xind=findgen(5001)
yind=exp(-xind/d0)
oplot,xind,yind,thick=8,color=19+doinstr
xyouts,2000,0.85,'Correlation decay length ='+string(d0,format='(I5)')+' km',$
  charsize=0.7
;
end