;
; Computes correlations between full, high and low pass timeseries of density
; instrumental anomalies.  Just does the ALL region.  Two instrumental
; timeseries are used: ALL and NHEMI.  Many ALL MXD timeseries are used, 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=1881.
peren=1960.
thalf=10.
;
; Get ALL instrumental timeseries and filter it
;
restore,filename='instradj_'+regname(ireg)+'.idlsave'
xall=x
yall=instradj
filter_cru,thalf,tsin=yall,tslow=ylall,tshigh=yhall,/nan
;
; Get NHEMI instrumental timeseries and filter it
;
restore,'instrts_nh.idlsave'
totval=total(nhts(*,3:8),2,/nan)
totn=total(finite(nhts(*,3:8)),2)
keeplist=where(totn gt 0,nkeep)
xhem=x
yhem=totval*!values.f_nan
yhem(keeplist)=totval(keeplist)/float(totn(keeplist)) 
filter_cru,thalf,tsin=yhem,tslow=ylhem,tshigh=yhhem,/nan
;
if total(abs(xhem-xall)) ne 0 then message,'Ooops'
;
; 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,'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))
;
; Open output file
;
openw,1,'corr_fixed.out'
printf,1,'Correlations between timeseries'
printf,1,'density vs. instru:         ALL               NHEMI'
printf,1,'Fixed-grid     Full  High   Low    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,$
  1720,1740,1760,1780,1800,1820,1840,1860,1880,1900,1970,1980,1990,1994]
ncut=n_elements(allcut)
;
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)
  ;
  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
  ;
  kl=where(finite(yall+ymxd) and (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) )
  ;
  kl=where(finite(yhem+ymxd) and (xall ge perst) and (xall le peren),nkeep)
  r4=correlate( yhem(kl) , ymxd(kl) )
  r5=correlate( yhhem(kl) , yhmxd(kl) )
  r6=correlate( ylhem(kl) , ylmxd(kl) )
  ;
  printf,1,ct,r1,r2,r3,r4,r5,r6,$
      format='(I4,7X,2(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