;
; Computes correlations between full, high and low pass timeseries of MXD
; anomalies.  Just does the ALL region.
;
; The idea is to see how the early, reduced sample, hemispheric mean captures
; variance.  An ALL timeseries is created from 1701 to 1960 using all
; available data (this is just the 1701 to 1960 section of the sample-size
; adjusted ALL MXD timeseries).  This is then compared to ALL timeseries
; generating based on fixed grids of those chronologies that go back to
; various times.  Since they are fixed grids and they are only needed as far
; as 1960, there is no need to apply a sample-size correction.
;
; Get regional information
;
restore,filename='reglists.idlsave'
ireg=0     ; 0=ALL
;
; Specify period and timescale separation to be used for correlations
;
perst=1701.
peren=1960.
thalf=10.
;
; Get ALL MXD timeseries and filter it
;
restore,filename='densadj_'+regname(ireg)+'.idlsave'
xall=x
yall=densadj
filter_cru,thalf,tsin=yall,tslow=ylall,tshigh=yhall,/nan
;
; Get all the MXD chronologies
;
ncid=ncdf_open('tree_dens_nh.nc')
ncdf_diminq,ncid,'time',dummy,ntime
ncdf_diminq,ncid,'station',dummy,nstat
ncdf_varget,ncid,'year',xmxd
ncdf_varget,ncid,'latitude',ylat
ncdf_varget,ncid,'longitude',xlon
ncdf_varget,ncid,'density',density
ncdf_attget,ncid,'density','missing_value',valmiss
ncdf_varget,ncid,'fraction',weight
ncdf_close,ncid
misslist=where(density eq valmiss,nmiss)
density(misslist)=!values.f_nan
if total(abs(xmxd-xall)) ne 0 then message,'Ooops'
;
; Extract those in the ALL region
;
mxdall=density(*,treelist(0:ntree(ireg)-1,ireg))
wtall=weight(*,treelist(0:ntree(ireg)-1,ireg))
xlon=xlon(treelist(0:ntree(ireg)-1,ireg))
ylat=ylat(treelist(0:ntree(ireg)-1,ireg))
;
; Open output file
;
openw,1,'corr_early.out'
printf,1,'Correlations between timeseries'
printf,1,'MXD vs. MXD:                ALL'
printf,1,'Fixed-grid     Full  High   Low'
;
; Now construct an MXD timeseries using all chronologies with data back
; to some cutoff time and filter it
;
allcut=[1400,1450,1500,1550,1600,1620,1640,1660,1680,1700,1800]
ncut=n_elements(allcut)
;
; Get ready for plots of distribution
;
multi_plot,nrow=4,ncol=3,layout='large'
if !d.name eq 'X' then window,ysize=850
loadct,39
map=def_map(/npolar) & map.limit(0)=30.
labels=def_labels(/off)
def_1color,cr,cg,cb,50,color='red'
def_1color,cr,cg,cb,51,color='blue'
;
for i = 0 , ncut-1 do begin
ct=allcut(i)
yw=where(xmxd eq ct)
yw=yw(0)
kl=where(finite(mxdall(yw,*)),nkeep)
print,ct,yw,nkeep
mxdsub=mxdall(*,kl)
wtsub=wtall(*,kl)
;
inter_boxfd,/nodata,map=map,labels=labels,title=string(ct,format='(I4)')
cpl_usersym,/circle,/fill
if i eq 0 then begin
  xsub=xlon(kl)
  ysub=ylat(kl)
  plots,xsub,ysub,psym=8,color=50,symsize=!p.charsize*0.6
  xold=xsub
  yold=ysub
endif else begin
  xold=xsub
  yold=ysub
  xsub=xlon(kl)
  ysub=ylat(kl)
  plots,xsub,ysub,psym=8,color=50,symsize=!p.charsize*0.6
  plots,xold,yold,psym=8,color=51,symsize=!p.charsize*0.6
endelse
  ;
  ; NOTE THAT WE CAN IGNORE MISSING DATA SINCE WE WANT A FIXED GRID!
  totwdens=total(mxdsub*wtsub,2)
  totwvals=total(wtsub,2)
;  totwdens=total(mxdsub*wtsub,2,/nan)
;  totwvals=total(float(finite(mxdsub))*wtsub,2)
  ymxd=totwdens/float(totwvals)             ; this is the weighted anomaly
  mknormal,ymxd,xmxd,refperiod=[1881,1960]
  filter_cru,thalf,tsin=ymxd,tslow=ylmxd,tshigh=yhmxd,/nan
  ;
; shouldn't have any missing!
;  kl=where(finite(yall+ymxd) and (xall ge perst) and (xall le peren),nkeep)
  kl=where((xall ge perst) and (xall le peren),nkeep)
  r1=correlate( yall(kl) , ymxd(kl) )
  r2=correlate( yhall(kl) , yhmxd(kl) )
  r3=correlate( ylall(kl) , ylmxd(kl) )
  ;
  printf,1,ct,r1,r2,r3,$
      format='(I4,7X,2X,3F6.2)'
  ;
endfor
  ;
printf,1,' '
printf,1,'Correlations carried out over the period ',perst,peren
printf,1,' '
printf,1,'To separate low and high frequency components, a gaussian weighted'
printf,1,'filter was used with a half-width (years) of ',thalf
;
close,1
;
end